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1.0 Introduction 

The goal of most researchers studying “how to see” is object recognition. The ap¬ 
proach includes separating a scene into objects, describing these objects in terms of 
some model such as generalized cones (Binford, 1971; Marr,l982), and then match¬ 
ing the result to a canonical description in memory (see Ballard k Brown, 1982). 
Surprisingly, the most progress has been made in the last two steps, whereas the 
initial first step of identifying objects in a scene is still largely unsolved in computer 
vision and is not understood at all for biological systems, except in special cases 
where an isolated object moves against a stationary background, or is distinguished 
in three-dimensions with no occlusions present. 

Our work differs from most in vogue today by suggesting that a scene be 
separated into objects not by edge extraction followed by contour grouping, but 
rather by recovering regions of differing material qualities. This approach was 
attempted early in scene analysis (Ballard k Brown,1982), but then was based only 
on simple image features such as color and texture. Although a useful beginning, 
a more reliable method is to use the image features to infer material properties of 
the world (e.g. water, cloud, tree, grass, wood, asphalt, etc.) and to break up the 
scene into its “stuff”. Such world-based properties are stable and survive changes 
in illumination, viewpoint, shading, surface orientation, etc., which confounded 
earlier image-based attempts at scene analysis (Heeger k Pentland 1987; Kube k 
Pentland, 1986). 

By identifying the “stuff” a region is made of, we also gain the additional 
benefit of knowing something about what object that region might correspond to. 
For example, hair implies an animal, feathers L'rd, grass a surface which is part 
of the terrain (Richards, 1982; Richards k Bod V 1988). 

Our major findings are : 

• Useful material descriptions must include both shape and reflectance predi¬ 
cates. (In many cases, dynamics is also needed.) Shape topology is far more 
important than precise metrics. 

• The Cook-Torrance reflectance model can be simplified to a three-dimensional 
system for purposes of rendition and image understanding. 

• The use of Graphics Psychophysics (i.e. image analysis by generating synthetic 
images) is a very powerful but time-expensive tool for image understanding. 

1 The frontispiece is a composite of computer-generated surfaces and pictures taken from 
Brodatz (1966). The textures for all these surfaces can be modeled by fractal processes. 
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a Angle between N and H 

6 Angle between L and H or V and H 

X Wavelength 

D Facet slope distribution function 
d Fraction of reflection that is diffuse 
du, Solid angle of a beam of incident light 
£, Energy of the incident light 
F Reflectance of a perfectly smooth surface 
/ Unblocked fraction of the hemisphere 
C Geometrical attenuation factor 

H Unit angular bisector of V and L 

/, Average intensity of the incident light 

/„ Intensity of the incident ambient light 

I, Intensity of the reflected light 

A, Intensity of the reflected ambient light 

k Extinction coefficient 

L Unit vector in the direction of a light 
m Root mean square slope of facets 

N Unit surface normal 

n Index of refraction 

£, Ambient reflectance 

R Total bidirectional reflectance 
Diffuse bidirectional reflectance 
R, Specular bidirectional reflectance 

j Fraction of reflectance that is specular 

V Unit vector in direction of the viewer 
w Relative weight of a facet slope 
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Figure 1 Parameters of the Cook ic Torrance (1982) model for surface re¬ 
flectance. 


• We present the first study in computational acoustics showing how material 
properties can be recovered from sounds. 

• More work is needed to understand the useful image correlates of the motion 
of non-rigid objects, especially those materials without surface markers such 
as rubber or water. 


2.0 The Problem 

Recovering material properties from image information is difficult because the avail¬ 
able signal, namely the image intensity, is a function of several factors which are 
confounded. For example, the observed intensity (flux) / of a point on the retina 
depends upon (l) the spectral reflectance (albedo) of the surface, (p(A)), (2) the 
strength and spectral composition of the illuminant E(A), (3) the geometry and 
structure of the three-dimensional surface N(<r, r) relative to the light source direc¬ 
tion L and the viewing position V, and finally (4) the reflecting properties of the 
material R. In particular, 
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I(x,y,\) = p(\)*E(\)*(NL)*R (l) 

Thus, if one wishes to “know” what kind of pigment p(A) makes up the sur¬ 
face (chlorophyll, xanthophyll, flavanoid, ferrite, etc.) the observer must somehow 
discount the illuminant E (including shadows), the surface orientation N, the illu- 
minant direction L, his particular viewing position V, and any specular component 
of the surface. (The latter component is embedded in R, whose complete descrip¬ 
tion requires many of the additional parameters defined in the table of Figure 1.) 
Obviously with so many unknowns, any single intensity value is insufficient (Horn, 
1977, 1987). Rather, several different regions in the image of an object must be 
examined for evidence that will disambiguate among all the possible origins of the 
image intensity values (Richards, 1988). For example, to tell whether an urn is 
made of brass, or rubber, one must “inspect it” and integrate information from the 
highlight, the occluding boundary at two or more locations, the intensities arising 
from shaded regions, etc. 

Considering the complexity of the interrelated unknowns, the problem of recov¬ 
ering the material property from image intensities seems impossible. Yet observers 
can indeed distinguish between materials such as metal, plastic, rubber, ceramic, 
etc., even when all other factors such as object shape and illumination are constant 
(see illustrations of Cook <k Torrance, 1982; Richards ii Ullman, 1987). On what 
basis are these inferences made ? 

3.0 The Approach: Graphics Psychophysics 

Obviously an observer can not recover all the physical parameters of a surface. Only 
some of these are relevant to the inference process. Our method for identifying these 
relevant parameters is simple: we use computer graphics to generate an image that 
“looks like” the real surface. 

Figure 2 illustrates in more detail the role “analysis-by-synthesis” can play 
in Image Understanding. When a scene is imaged upon our retina, it becomes a 
pattern of receptor activities, as shown in the lower right of the figure. This is what 
we are given, namely a 2D array of numbers. Our task is to infer the structure and 
stuff of the world, as shown by the Inference arrow on the upper right. We face 
two problems. First, we need to know the numbers in the array. Second, we need 
to know how the numbers “got there”. 

For any given scene, finding the numbers in the image array is trivial. Simply 
take a picture with a vidicon and read the numbers out of frame buffer. We can 
check to be sure these numbers are correct by displaying the contents of the buffer 
to see if the 2D display looks like the reed 3D scene. In terms of Figure 2, the 
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Figure 2 The role of computer graphics. 

pattern of numbers in the display’s buffer then matches that of the image of the 
scene (or at least there is a one-to-one mapping). 

The tough problem is the second one, determining why the 3D surfaces in the 
world led to the particular array of numbers. Without knowledge of this mapping, 
the inference problem can’t be solved, for we have not specified the function we 
wish to invert. In terms of our figure, we need to know how to load the frame 
buffer to duplicate the real 3D scene. This is a problem in computer graphics. 
The generation of a realistic-looking cloud or mountain forces the programmer 
to model the world, to make explicit the optical properties of a surface that are 
psychophysically relevant (Blinn, 1977, 1982). The more realistic the display, the 
closer the programmer’s model captures the relevant properties and geometery of 
the surface. Specifically, we know the origin or cause of each image intensity element 
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/(*, y, A, t) on the retina of our eye, for we have specified function explicitly in our 
program. Hence the mapping from the world to the image plane is known. It then 
becomes apparent what external physical properties underly the observed property 
of “metal” or “rubber”, as in the Cook &i Torrance reflectance model (1982), or 
the rough surface of a mountain, as in the fractal models of Carpenter (1980), Voss 
(1981) and especially Pentland (1984). Computer graphics thus provides us with 
the opportunity to complete our understanding of the mapping of the world into 
the image. Once understood, possible inverse mappings and constraints can be 
intelligently explored. Inferences from images can proceed. 
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4.0 Inferring Material Types 

In the sections to follow, we illustrate our Graphics Psychophysics approach with 
two vision examples, inferring “wood” and “water”. Appendix II also includes 
an analogous study of how a natural property might be recovered using acoustic 
information. 


4.1 Why Wood Looks Like Wood (D. Honig) 

1.0 Introduction 

The texture and pattern of a surface provides a major source of visual information 
about the material or “stuff” we are looking at. Once the type of material is 
classified, a host of other facts can be accessed, such as the expected rigidity of the 
surface, its hardness, its density—or other properties that may guide us in object 
recognition or manipulation. For the human observer, the classification of materials 
using optical information is almost automatic and consequently seems trivial. Do 
we simply store a vast collection of textures typical of surfaces of interest, learning 
by example the proper associations? Or, alternatively, do we first infer from the 
optical array the type of process which created the pattern, and then proceed with 
our inference? Here, we take the latter approach: how can the process which 
created a visible surface structure be inferred? We illustrate this approach using 
the pattern created by wood grain. 

2.0 The Method 

To understand why wood looks like wood, we begin by creating a rendition of wood, 
using computer graphics. This rendition is not just a texture mapping or “painting” 
of a pattern onto a surface, but is actually a three-dimensional construction based 
upon tree structure and how trees are cut up into lumber. lYees are a natural 
category of plant with certain family traits. Common structured properties of trees 
combine with the ways in which lumber is commonly cut to yield a pattern of 
textural elements that is practically unique to wood. It is these kinds of regularities 
that allow us to classify things. By using computer graphics to create a 3D model 
of the tree structure and then by cutting this model structure, we can identify the 
origin of each regularity typical of the wood grain pattern. 
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Figure 1 The structure of a tree (from Sligo, 1985). 


S.O The Proeeaaea 

As previously mentioned, the pattern typical of wood grain depends both upon the 
natural structure of trees and upon the way trees are cut into lumber. We discuss 
each process separately, prior to presenting our model of these processes. 

S.l The Tree 

A tree is a living process, absorbing minerals from the ground and photosynthetic 
energy from the sun. The trunk and its branches are cylindrical structures with 
three main parts: the cambium, the pores, and radial compartments (see Figure 1). 

S.l.2 The Cambium 

The cambium is the only live part of the tree aside from the leaves. It is a thin, 
vascular layer in the trunk and in all branches and roots. Bark protects the cambium 
on the outside. (Here we neglect the bark, which is stripped prior to cutting 
boards.) The cylindrical cambium shell grows outward by depositing concentric 
layers of xylem (wood) cells on its inner surface. Fine tubes that transport water 
and other minerals up the tree are interspersed throughout the woody tissue. On 
the outer surface of the cambium layer the phloem (inner bark) is deposited, which 
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Flgnre 2 Typical ways in which trees are cut into lumber. 


transports photosynthetic products down the tree. The outer surface of the inner 
bark is the bark which protects the tree. 

The annual “growth rings” of a tree are formed because as the cambium grows 
outward, depositing new wood on its inner surface, the new tissue gets darker 
toward the end of the growing season. As winter approaches, the tree accretes new 
wood more slowly, so the natural pigmentation is denser. 

8.1.8 The Pores 

The tree transports various fluids through pores. These pores are simply cellulose- 
lined cells or fine, tubular structures oriented parallel to the axis of the trunk or 
branch. 


8.1.4 Radial Compartments (Rays) 

Trees defend themselves against infection by partitioning, walling-off damaged and 
infected tissue. Two of the interior walls of the partitions are radial planes of 
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Figure 3 Several special cuts of a tree into boards. The shrinkage charac¬ 
teristics of lumber depend on how it is cut (from the U.S. Forest Products 
Laboratory). 


parenchymal tissue distributed among the woody tissue. Vertically, trees stop the 
spread of disease by blocking the tubes (pores) passing through the site of injury. 
Along the radius, trees wall-off infection with the annual growth rings, whose high 
density serves as an effective barrier. The compartments are most visible in cross 
section, for wood will often crack along the radial walls of compartments. 


3.2 Lumbering 

When trees are cut into boards, there are typically only three major types of cuts: 
the “stump” cut; the “board” cut; and the “plywood” cut (see Figure 2). 

Let the axis of the trunk or branch be oriented in the y-direction. A “stump” 
cut is simply a plane whose normal vector is parallel to Y. The typical “board” cut 
is a plane parallel to Y, i.e. parallel to the axis of the branch. These cuts can have 
rather complicated cross-sections, as illustrated in Figure 3. The “plyboard” cut is 
a cylindrical shell or “stripping” parallel to the cylindrical structure of the trunk or 
branch. Finally, there are several “odd cuts”, like spheres or dowels. Most of our 
renditions used board cuts, which generate the most common wood grain patterns. 
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Figure 4 Simplified model of the three major components of a tree’s structure. 


4.0 The Model 

Our idealized tree has three structures: a series of cylindrical shells (cambium), fine 
tubular vessels running parallel to the cylindrical axis (pores) and radial compart¬ 
ments (rays). As illustrated in Figure 4, these we identify as the “rings”, “rays”, 
and “pores”. 

The “rings” in the computer model were created by concentric cylinders with a 
fixed shell width, but having a gray scale gradient across the shell. The pores were 
simulated by randomly distributing tiny vertical tubes throughout the modelled 
tree trunk (i.e. the nested cylinders). The rays were modelled as radial vertical 
rectangular planes of bounded radial and vertical extent. These planes were ren¬ 
dered darker than the backgound material, spaced at roughly 20° intervals. 

Although wood is usually light brown, some varieties like “rosewood”, “cherry”, 
or basswood can exhibit a range of pastel colors. Furthermore, if wood is stained 
or n nted, it can be of almost any color. Here, we generally used 256 brightness 
le-/t • of the same color. One part red to one-third green to one fifth blue in the 
C rmonitors produced a realistic wood-brown. 
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5.0 The Renditions (Instant Psychophysics) 

Each rendition can be characterized by the plane of section through the tree trunk 
(or branch). With the K-axis of the coordinate system aligned with the central axis 
of the tree, and with the Z-axis pointing to the viewer, any plane can be described 
by its normal vector [x, y, z]. 

Figure 5A shows a “stump cut”, N = (.? 1 .1], overlaid with the tapered end of 
a dowel, N — [—110]. For the stump cut, we see the rings manifest as circles, with 
the rays appearing as radial lines and the pores as points. We have also added a 
slow radial modification in albedo to model longer-term variations in density due 
to climate. This addition breaks up the otherwise regular geometric pattern and 
produces a more acceptable rendering. As an aside, we note here that the particular 
waveform of the annual density of wood is not very critical (i.e. the cross-sectional 
intensity profile of the rings). Technically the most dense (and hence the darkest) 
portion of the ring should be on the inner surface of each layer of the cambium. In 
practice, the typical observer ignores this subtlety—we have reversed the phase of 
the pattern with no loss in rendition quality. 
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The remaining figures show “board cuts” N = [a 0 6], with a, b taking on 
slightly different values for each of the five panels of each figure. The most ob¬ 
vious feature of all these panels is the expected elliptical pattern of the rings. In 
Figure IB, only the rings (plus slow climatic modulation) is illustrated. The ren¬ 
dition is acceptable. Adding the rays greatly improves the rendition (Figure 1C). 
However, breaking up the regular elliptical pattern by modulating the cut produces 
the most realistic wood grain. Note that an undulating plane of section is equiv¬ 
alent to a planar section of a tree which is not strictly a nested stack of regular 
cylinders. The observation that the presence of irregularities is a necessary part of 
a good rendition is important. It means that the inference “wood” requires some 
evidence that a strictly regular, machine-like process is not solely responsible for 
the observed structure. 

In sum, our renditions suggest that the geometrical model consisting of nested 
cylinders, radial compartments, and pores (which create a fine-grained “peppery” 
texture) is acceptable provided that two kinds of perturbations are included in the 
strictly geometric model: 1) longer-term climatic variations in pigmentation density 
and 2) allowance for undulations in the cylindrical axes and diameter. 

6.0 The Inference Process 

A wooden surface typically exhibits a “grain direction”, either from the rays or from 
the pores or both. For perfect planar cuts, this grain direction will range from radial 
(stump cut) to parallel (board cut). The major axes of the ellipses of the growth 
rings must align with this grain direction, all intersecting together at the center of 
the tree. In addition, the spacing between the elliptical rings must be the same 
along any given ray (or grain direction). Furthermore, the rings must exhibit an 
albedo variation whose function is normalized for the different separations between 
the rings for different rays. This is the strictly geometric model. In addition, some 
evidence must be presented that the pattern is not simply a regular geometric 
structure. Slow albedo modulation and perturbations in the shape of the ellipses 
along a given ray seem to provide convincing evidence of a natural growth process. 1 


7.0 False Targets 

One simple approach to enumerating the important false targets is to recognize that 
the strictly geometric model results in an image of concentric ellipses with equal 
radial increments from one to the next ellipse. We now have four free parameters 
as this plane is moved to create a nesting of surfaces in 3-space: the ratio of major 

1 Note: saw cuts may also create density variations parallel to the edge of the boards. 
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to minor axis, the length of the major axes, and two degrees of freedom in the 
direction of movement of the common axis (i.e. center of the ellipsis). 

An onion falls into this class of possible generating surfaces. In this case, the 
major and minor axes of the ellipse are equal (circles), and the plane of nested 
circles is moved in a direction perpendicular to the original plane of the circles, and 
the major axis (radius) is decremented as we move away from the initializing plane, 
yielding a set of nested circles (or ellipses if we wish). A nesting of cones also could 
yield a false target. 

To exclude these possibilities, we may require further evidence, such as radial 
“cracks” or albedo variations in the ring pattern. But the key reason for the success 
of the inference “wood” from the image slice is that in the natural world, structures 
created from nested cones or ellipses are very, very rare—indeed it is difficult to 
give a natural example of such a structure! A nesting of spheres is possible, such as 
the onion, but here the potential for confusion is only with the less common “stump 
cut”. Fortunately, Nature creates structures according to rules which tend to group 
similar classes together, creating “modes” of particular structures rather than a 
continuum of all possible functions (Bobick, 1987). The cylindrical construction of 
plants and trees is an example of such a “mode” which permits successful inferences 
from the pattern we recognize as “wood”. 


8.0 Summary 

One might ask, was computer graphics necessary to determine the characteristic 
features of wood? Couldn’t we have just sketched the projections of the main 
components of the tree projected onto the image plane, just as an artist might do? 

Certainly an artist’s rendition of any scene is informative and helpful. But 
can the artist tell you exactly what each brush stroke represents in terms of a 
world property? Without such knowledge, the inference from image to world is not 
possible. Furthermore, the artist is depicting his internal model of the material, 
which is not necessarily an accurate world model—certainly not a physical model. 
Often, certain features are exaggerated for emphasis, and usually only common 
views are depicted. The graphics model, on the other hand, is not limited to special 
views and exaggeration must be made explicit. Hence the graphics model can be 
used to explore unusual viewing or illumination conditions to determine the actual 
character of the observer’s internal model. Such power is needed if the observer’s 
representation and inference process is to be understood fully. Further examples 
of this point follow in subsequent sections, where unusual viewing conditions of 
a “correct” physical model lead to mis-perceptions because the observer’s model 
assumes particular viewing and illumination conventions. 
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4.2 Inferring “Water” from Images (T.J. Kung) 

1.0 Introduction 

Water waves are recognizable from a black and white photo. What is the image 
information that allows us to make this inference? In order to answer this question, 
we use computer graphics to simulate water. Once a perceptually deceiving sim¬ 
ulation of water is obtained, then the study of the inference process cam proceed. 
This report presents a model for generating near-natural images of water waves. 

Computer generated water waves have been created previously for flight sim¬ 
ulators. However, these images did not look very realistic, probably because very 
simple models were used. For instance, Max (1981) presented a ray-tracing pro¬ 
cedure in which ocean waves and islands are rendered. Although he was able to 
render good reflection images off waves, the ocean waves do not look like real waves. 
Ogden (1985) generated a real-time animation sequence of water by using the Burt 
Pyramid (1983), but again her model was not based on the physics of water. Fur¬ 
thermore, because the method essentially generated fractals in the image plane, 
there was no real description of the 3D surface. Here, our objective is to create a 
3D surface of water undergoing wave motion, and to be able to view this surface 
from any arbitrary angle. Complicated effects like the reflection of an object off the 
water surface are not considered in this paper. 

We divide the problem into three parts: (1) reflection, (2) surface shape, and 
(3) dynamics. In interpreting images, basically we are trying to infer the surface 
properties and shape from the image intensities. As a first step in understanding 
this inference process, we must understand what properties make water charac¬ 
teristically different from other surfaces. If we knew these properties and how to 
describe them mathematically we should be able to create a realistic image. Among 
the surface properties, we are only interested in the optically related properties, i.e. 
the reflectance function. A dynamic view yields a considerably stronger impression, 
therefore a brief discussion of wave dynamics is included also. 


2.0 Reflectance Function 

Because water is a transparent material, the intensity and color of the visible light 
emerging from the sea depends on 

(i) The illumination. The elevation of the sun and the sky and cloud cover deter¬ 
mine the quantity and spectrum of downwelling incident light. 
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Figure 1 Synthetic water images are shown in the two left-hand panels. The up¬ 
per right panel is an artist’s rendition. A photographic image of water is shown 
in the lower right. Most people prefer one of the synthetic images, demonstrating 
that our model is capturing the psychophysically relevant information. 


(ii) The optical properties of the sea-water itself. The light which penetrates into 
the sea is spectrally modified by absorption inside the medium before being 
partly backscattered toward the atmosphere. These absorption and backscat- 
tering processes have their origin in the water molecules, and also in all the 
other substances present in the sea as dissolved or particulate matter. 

People can recognize water waves from a black and white photo, therefore the 
most important thing in perceiving water waves is the relative intensity of emerging 
light from different places of the water surface, not the details of spectrum and 
quantity. This fact makes it possible to simplify the reflectance function without 
considering the functional dependency on the wavelength of lights. 

Cook k Torrance (1982) proposed a reflectance model for computer graphics 
as follows: 

I — RambEamb + ^ ^ •®dir(N • L)((ff2 ma i + sR t pec) (l) 

where E is the intensity of light source, N is surface normal unit vector, and L 
is the light direction of a point light source. The point light sources are linearly 
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Figure 2 The geometry of reflection. 

combined with the ambient light source. Each point light source is further divided 
into a specular component and matte component. The d and a in equation (1) are 
the fractions of reflectance that are matte and specular. The specular component 
is given by: 

Rtpec = 7 (N • L)(L • V) (2) 

where Fisa Fresnel term which describes how light is reflected from each smooth 
microfacet. It is a function of incidence angle and wave length. The geometrical 
attenuation factor G accounts for the shadowing and masking of one facet by sm¬ 
other. The facet slope distribution function D represents the fraction of the facets 
that are in the direction H (Figure 2). 

In rendering water, we use Cook k Torrance’s model with some modifications. 
As discussed before, water backscatters part of the incident light. Although the 
backscattered light is a function of the incident light spectrum, optical properties 
of the water itself, and the elevation of viewing position, we assume these condi¬ 
tions are fixed in a limited area and within a short period of time. Therefore the 
backscattered light is a constant and is uniform over the water surface. 

Figure 3 shows the value of Fresnel term for various indexes of refraction n. 
The Fresnel term is significant only when the incident angle is greater than 70 
degrees. For water the index of refraction equals 1.32. Based on this property, we 
make two simplifications. First, we neglect the matte term because a water surface 
only reflects light in very limited directions. Second, we neglect the point light 
source-sun. By doing so we eliminate the special cases when it is dawn and when 
it is dusk, and also scintillation effects which are easy to mimic. 

The ambient light source (sky) is considered as a hemi-sphere light source. 
No matter where the viewing position is, it is always possible to see some lights 
from somewhere in the dome directly reflected off the water to the eyes. When 
the viewing angle is less than 70 degrees, so is the incident angle, and the reflected 
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Figure S Fresnel reflectance as function of incident angle. The parameter on 
the curves is the index of refraction, n. Note the change in scale for the right 
portion of the graph. 


light is much less than when the viewing angle is greater than 70 degrees. In other 
words, the intensity of reflected light from the water surface is a function solely of 
the angle between the viewing direction and the surface normal of the facet. 

As a conclusion of above discussion, our reflectance model for water waves is 

Aotai = 4. + 4m6 J^' (N-L)(L-V)^ ^ 

where 4« is the backscattered light intensity and is a constant. The integration 
term is for the ambient light source. Since this term is a function of the angle 
between N and V only, we calculate these values numerically for different angles 
and store them in a table. In rendering the images, we then just find the value of 
this integration term from that table. In order to simplify the integral work, the 
Gaussian model is used for the facet slope distribution function D 

D = exp {-(a/m) 2 ) (4) 

where a is the angle between N and H (Figure 2). 0.8 is used for m. With this 
value D is small when a is greater than 20 degrees and then the integration can be 
neglected. In other words, the integrated region can be reduced from a hemi-sphere 
to a cone of 20 degr e radius. The expression for G is 
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~ ... . 2(N - H)(N • V) 2(N • H)(N • L) 

-(V~H)- 1 -(V"Hj- 


(5) 


The cone is further divided into three areas based on which term in Eq. (5) is 
important. The full expression of F used in the integration term is 


1 {g-c? f \c(i + c)- i] 2 ' 

2 (f + c) 2 |c(y + c) + l] 2 . 


( 6 ) 


c=V H 
, 2 =« 2 + c 2 -l. 


For water, n is 1.32. The ratio of ht/Iamb was set to 1/3 in our simulation. 


S.O Surface Shape 

Traditionally, sinusoidal functions are used to model ocean waves (Bigelow k Ed¬ 
mondson, 1952; Gross, 1972). There is a strong theoretical basis to assume that 
the ocean surface is formed by sums of long crested sinusoids: 

= + + (7) 

Since this model requires the sum of a large number of wave forms to obtain a 
realistic picture, it is not practical for image generation. Schachter (1980) proposed 
a narrow-band noise model to overcome this problem. In this model, it is assumed 
that envelopes and phase shifts are slowly varying: 

w(x) = a{z) * cos (nz + ^(z)) (8) 

Although this model has the advantage of being able to simulate different textures, 
the results show that the water generated by this model does not have the shape 
randomness that real water does. 

We propose another model which is a summation of squared of sine waves with 
slowly varying wave length: 

w(*,v) = £[-A* ^1 + sin ( A ^^ (xc°s/?-t-ysi n i g) -t-^))] (9) 

Here we assume waves traveling in the x, y plane and 0 is the angle between x-axis 
and the traveling wave direction. The ideal wave form is sinusoidal but a lot of 
factors earn distort waves from their ideal forms, for instance winds, nonlinearity, 
and nonuniform water depth for shallow water waves (Bascom, 1980). Normal wave 
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Component Number 

1 

2 

3 

4 

5 

Mean wave-length l 

25 

30 


10.7 

16.7 

Amplitude A 

0.3 

0.4 


0.15 

0.25 

Modulation Rate r 

Bums 


IQilEH 

0.0 

0.0 

Angle (dgr) 0 

110 

90 

65 

120 

60 


Table 1 Static wave parameters. 


crests tend to go steep before collapse. Consequently, the distorted wave form is 
steep at the crest and fiat at trough. We found it takes less computation to get the 
distorted wave form by using a function which has its shape close to the distorted 
wave form such as the sine squares than by using Fourier series (summation of 
sinusoidal waves). [See Appendix II.] 

Because of nonlinear effects, wave trains are always modulated when traveling 
in & field. Instead of assuming slowly varying amplitude and phase shift as Schachter 
does, we assume slowly varying wave length. The amplitude is assumed to be a 
constant because the variation is too small to be noticed in short distances. 

Nonlinear instability theory tells us that instability causes steep gravity waves 
to evolve into the three-dimensional spilling breaking waves. This three dimensional 
wave form is the normal wave form we Bee in an open sea. To simulate the three 
dimensional wave patterns, it is necessary to have at least three components in 
equation (9). It was found that to generate a nice looking water wave image, each 
component should travel in a direction about 20 to 30 degrees apart from the other. 
The image looks even better with one more minor component on each side of the 
three major ones (with amplitude and wavelength about one third of the average 
values of the major components). 1 In Table 1 these directions are specified by 
parameter 0 which is the angle measured counterclockwise from the horizontal axis 
(x-axis). After testing several functions, we decided to use the following function 
for the slowing varying wave length 


A(z,y) = 


t 

1 + rsin(px + j>y) 


( 10 ) 


where / is the mean wave length and r is the normalized deviation of wavenumber, 
and p is the parameter which shows how slowly the wavelength varies. In our 


1 Texture psychophysics has shown that only 5 or 6 components are required to mimic 
very complex sums of sinusoids (Richards, 1979). 
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simulation p was fixed at 0.1. The actual values of the other parameters are listed 
in Table 1. The unit for wave length and amplitude is the pixel. 


4-0 Wave Dynamics 

The wave dynamics are solved from Navier-Stokes equation (Currie, 1975). With 
the assumptions that the Coriolis force are negligible, water is inviscid and incom¬ 
pressible, density is constant, and water motion is irrotational, the Navier-Stokes 
equation can be reduced to 

V 2 i> = 0 ( 11 ) 

where ^ is the velocity potential. There are two boundary conditions at the air- 
water interface. The kinematic boundary condition says that particles at free sur¬ 
faces remain on the free surface and the dynamic boundary condition says that the 
forces balance at free surface. One more boundary condition says that the vertical 
component of Quid velocity is zero on the ocean floor or bottom. 

The surface tension can not be neglected from the dynamic boundary condition 
when the wave length is less than 7 cm. These kind of waves are just the ripples on 
top of the larger gravity waves. By gravity wave, we mean that the gravity force is 
the only restoring force for a wave system. 

For a gravity wave, if the wave length A satisfies 

\ > 0.28A (12) 

where h is the ocean depth, then this wave is called a short wave or deep water 
wave. In this case, a wave travels with velocity 

C = (13) 

Because this velocity depends on wave length, longer waves travel faster than 
shorter waves and the wave system is dispersive. On the other hand, if the wave 
length satisfies 

k < 0.07A (14) 

This wave is called a long-wave or shallow water wave. In this case, the waves travel 
with the same speed 

C-v/F* (15) 
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Component Number 

1 

2 

3 

4 

5 

Mean wavelength l 

25 

30 


10.7 

16.7 

Deep water (h > 0.28A) 

mu 

MEM 



0.74* 

Shallow water (h < 0.07A) 

1.0* 

1.0* 

1.0* 

1.0* 

1.0* 


Table 2 Phase shift * between frames used to generate dynamic water waves. 


Therefore waves are nondispersive in shallow water. 

In order to have an idea about how these two modes of waves look, we simulated 
traveling waves in both cases based on our static model described before. We 
generated 40 static pictures, being the maximum allowed by our computer memory. 
Each picture was phase shifted from the previous one. For deep water waves this 
phase shift was proportional to the square root of wavelength and for shallow water 
waves it was independent of wavelength. The exact values of phase shift we used are 
listed in Table 2. We cycled these pictures sequentially and video taped the whole 
procedure. It turned out that the deep water mode looks better than shallow water 
mode. This result is consistent with our daily life experience. Except for shore 
waves and turbulent streams, we hardly see a pond, river, or sea with water depth 
less than a quarter of the average wave length. 


5.0 Conclusions 

A procedure for generating near-natural water-wave images by a computer has been 
presented. Three factors play the dominant role in our perception of static water 
images. The first is the Fresnel term in the reflectance function. The second is 
hemispheric illumination. Although a point light source, such as the sun, was not 
included in our simulation, this is obviously the key factor for the speckle and 
scintillation effects of water (Longuet-Higgins, 1960). The shape and composition 
of the wave-itself was also found to be important. For dynamic water-waves, it is 
necessary to include non-linear terms that mimic the behavior of deep water wave 
systems. 


6.0 Appendix I: Image Features 

In the main body of this paper we used computer graphics to determine which 
physical properties were significant factors in creating images that looked like water. 
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Fignre 4 The top panel shows a cross-section through a surface formed by 
two orthogonal sinusoids of equal wavelength and amplitude. The cutting plane 
contains the line of sight and two vertical surface normals, and hence is a special 
view. The cross-hatched region will appear as a darker “blob*, defined by an 
isoluminance contour corresponding to an emittance angle of 60°. In the lower 
panel, the surface normals are mapped on the Gaussian sphere, with the viewer 
direction normal to the plane of the page. The line of aero Gaussian curvature 
is shown as dashed (* = 0). 


Ignoring dynamics, these were i) Fresnel reflectance, ii) hemispheric illumination, 
and iii) a wave packet composed of (non-linear) sinusoids of low amplitude. Here 
we proceed with the second step of our graphics approach to image understanding, 
namely identifying an image feature characteristic of “water” or other fluids in 
motion. 

The top half of Figure 4 shows a cross-section of a water wave viewed from an 
angle of about 20°. The visual ray is tangent to the wave at R, which is a point 
on the “rim” of the surface (Koenderink ii van Doom, 1982). The locus of all such 
rim points will be a great circle on the Gauss map. If we now place the viewer 
direction normal to the page (bull’s-eye), then the rim will be the outline of the 
Gaussian sphere, as illustrated in the lower half of the same figure. 

Consider now the isoluminance contour of 5% reflectance on the wave surface. 
For a Fresnel reflectance function, this corresponds to a 60° viewing angle. We 
wish to determine the general shape of this contour, and hence of the dark region 
it encloses. Referring to the Gauss map, the region of 60° Fresnel reflection will be 
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described simply by a circle centered about the bull’s-eye, assuming hemispheric 
illumination without wavelet occlusions. 

To determine which part of the 60° Fresnel circle is visible, consider the locus 
on the sinusoidal surface which has zero Gaussian curvature (/c = 0). For our 
surface waveform of crossed sinusoids, this locus can be shown to be a distorted 
square (or rectangle) on the Gauss map, with the square centered about N. The 
region contained by both the Fresnel circle and the k = 0 square thus represents 
the surface orientations of the dark “blob". On the Gauss map, this shape will have 
cusps at the extremal views. (Obviously the exact choice of cut-off for the angle of 
incidence is not critical.) 

If the “blob” region on the Gaussian sphere is now mapped back onto the 
surface, we obtain two dark regions abutting at the locus k = 0. We wish to 
show that the extremal ends of this region, namely where the surface normal has 
orientation M, will be seen as “cusped”. (A cusp on the Gauss map in itself is not a 
sufficient condition for a cusp on the surface.) One approach is to solve analytically 
for the isoluminance contour and show that the first derivatives are discontinuous 
here. Alternatively, we note that as we move toward M along the blob contour, for 
there not to be a cusp at M , the isoluminance contour must be perpendicular to the 
parabolic line k = 0, because M is an extremal point of the contour. But generally 
this will not be the case, because the behavior of the isoluminance contour and the 
K = 0 contour are controlled by independent processes, and hence generally will 
intersect transversally as they do on the Gauss map. The extremal end-points of 
the dark “blobs” of water should thus be cusped. (Close inspection of photos of 
water confirm this prediction.) 

In the above, our analysis assumed a special viewing position and a special 
waveform. First, it should be obvious that the cusped property of the blob will 
remain unchanged if the wavelength of either one of the sine waves is increased 
(except in the limit as the wavelength becomes infinite). Similarly, we expect no 
effect of small non-linearities of wave shape. However, it is clear that viewing angle 
can alter the geometry of the blob drastically, for example when the surface normal 
N approximates the viewer direction, such as when we look almost directly down 
at water. In this case the tangent to k = 0 at M may lie in the plane of the view 
vector, causing the cusp angle to reach x/2, thereby eliminating the discontinuity 
at M. The cusped feature for water is thus guaranteed only for shallow (< 30°) 
viewing angles. 3 

One final point can be made about the shape of these dark “blobs” character¬ 
istic of water. Once the amplitude of a wave is such that its bitangents intersect at 

3 Note that for very shallow view angles the lower portion of the dark blob of one wave 
will be occluded by the bright ridge of the next nearer wave, also resulting in cusped 
ends, or even in an extra pair of cusps, depending upon the view angle. 
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Figure 5 A “zoom” photo of a water scene, showing clearly the predicted cusps 
and image features (credit for photograph - unknown). 


less than 120° (see Figure 4), then the wave will crest. Because the wave silhouette 
(rim) sets bounds on the shape of the blob, it too must be elongated sufficiently in 
order not to violate this constraint. Our characteristic image feature for water will 
thus be a very elongated, horizontal elliptical-like blob having cusped end-points at 
its extremal viewing positions. It is of interest to note that Ogden (1985) created 
an illusion of a water surface by convolving with image noise a mask consisting of 
two such blobs of opposite sign (simply by splitting an ellipse). Similarly, properly 
chosen packets of Gabor filters can also generate “water-like” textures by (thresh- 
olded) convolution with noise (Daugman, 1988). Both of these masks create image 
patterns with dark, elongated, cusped blobs. 


7.0 Appendix II: Wave Shape 

Here we illustrate the psychophysical effects of three different choices for the shape 
of water waves. 

In the top row of Figure 6 the non-linearity in the shape of each wave is created 
simply by squaring a sinusoid (this was our final choice). 
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Figure 6 The psychophysical effect of the shape of the water wave upon its 
appearance, using Fresnel reflection and hemispheric illumination. 


For comparison, the next row uses a square of a triangular wave plus its second 
harmonic. Note the negligble difference in the appearance of the two types of waves 
under fresnel reflection and hemispheric illumination (top and middle right). 

Finally, the last row shows the appearance of water if only the linear sinusoid 
or the fundamental component of the triangular waveform are used. Neither of 
these renditions looks as good as the non-linear choices above. 
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5.0 A New Reflectivity Model 

Spectral reflectance, or more properly albedo p{ A), is probably the simplest, and 
certainly the most colorful surface parameter in the image-intensity equation. But 
in addition, we must also consider R —the reflectance function itself. 

I = p{\) E(\) (N L) R (1). 



Figure 1 The underlying surface is two orthogonal sinusiods. The illumination 
and reflectance conditions can drastically alter the inferred shape and material 
type. Note that special conditions are required to create an image that looks 
like water (lower right). 
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Equation (1) summarizes again our basic relation between surface parameters, 
p{ A), N ; the illumination E{\), L, and image intensities. Cook &c Torrance have 
elaborated this equation in detail, with special emphasis upon the reflectance pa¬ 
rameter, R. Here, the matte versus specular aspects of the surface are collected, as 
well as textural effects. Following Torrance ic Sparrow (1967), the textural com¬ 
ponent of the surface is considered only at the optical level, so the micro-texture 
itself is not visible [although see Blinn (1978) for a treatment of visible textural ef¬ 
fects]. Although the interaction between the illumination conditions and reflectance 
parameters is obvious from equation (1), the full force of this interaction did not 
strike us until we began to model water. Furthermore, as shown in Figure 1, the 
visible surface structure, not just the invisible microfacet distribution, is critical to 
how we infer surface material type. In this figure, the underlying surface structure 
for all six panels in the top two rows is simply two orthogonal sinusoids at 90° to 
each other. The upper left panel, which uses Phong shading with some diffuse (i.e. 
matte) illumination, looks roughly “correct”, as does the Cook <Sc Torrance rendi¬ 
tion using an index of refraction of 10. Note the importance of adding the matte 
component to a surface, for when this is taken away (upper right), the underlying 
surface can not be seen. Note also that the choice of index of refraction is not crit¬ 
ical as long as there is a dominant matte (diffuse) reflectance component (compare 
left-most in second row with top-middle in first row). The refractive index does, 
of course, critically affect the specular reflectance (middle, second row versus top 
right, first row). Finally, rather than having just a single light source illuminating 
a completely specular surface (top-right or middle, second row), let the source oc¬ 
cupy the entire hemifield (right-most in second row). The result looks more like 
ribbon candy or an E-M photo than crossed-sinusoids. Yet once the amplitude of 
the sinusoids is reduced and made typical of water, the surface indeed begins to 
look like a common one—namely water (lower right). Clearly there is an important 
interaction between the visible structure of the surface, its reflectance properties, 
and the illumination typical for that surface. 

To capture this perceptual effect, we have recently proposed a new, psycho- 
physically-based reflectance model that makes the visible surface roughness an ex¬ 
plicit parameter of the representation of material type (Richards, 1988). 

The basic idea is that facets on the surface may appear over a range of scales, 
some visible, some not, depending upon the viewing distance. What is needed is 
a parameter that captures this surface property. One natural choice is a fractal 
measure. The new representation includes the parameter h (roughness), and q 
(specular component), as illustrated in Figure 2. Each takes on values between 0 
and 1. For example, whenever p, h or r\ is zero, then the surface will appear black, 
for then (i) either the albedo p is 0, or (ii) the surface is so rough that no light can 
escape to the viewer (a black body!), or (iii) the surface is a perfect mirror with 
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Figure 2 Representation for reflectance. Axes are: albedo, p\ Hansdorf dimen¬ 
sion less 1, h, and a measure of the index of refraction, r), scaled to be between 
0 and 1. 


index of refraction t} = 0, where light can not be reflected to the viewer except 
under very special conditions (such as retro-reflect materials). Intuitively, surfaces 
with a high fractal dimension (large Hausdorf-Besicovitch number) will be “rough” 
and hence will be dark regardless of albedo (p) and refractive index (q). Whereas 
surfaces with a low fractal dimension (H = 1 - K) will be smooth, and can be 
either dark or light depending upon (p) or (rj). This insight gives new meaning 
to “facet distribution” for both matte (p) and specular (»j) reflectance. As shown 
by Pentland (1983, 1984, 1988) this measure can be recovered easily from image 
information for matte surfaces. Appendix I extends his proof to glossy surfaces. 
Hence a fractal facet distribution function can be estimated from images, at least, 
for surfaces which approximate fractals—a very large class. 2 

Thus, in our formulation for the reflectance function, all the matte properties 
of the surface are captured by the albedo parameter, p. All the specular properties 
thus appear in the reflectance function, R, whose two major parameters are rj (in¬ 
dex of refraction) and h, which describe the surface roughness. 3 The recovery of 
the surface material would require estimating these three parameters from image 
intensities. As already noted, this is a difficult task in the general case because of 
the complexity of the full image intensity equation (1). However, in many cases, 

2 As an alternative to the fractal measure, see Lewis, 1987. 

3 Note also that X may have an orientation component, such as water does (i.e. approxi¬ 
mately 1-D) whereas clouds are typical of a 2-D fractal function (see frontispiece). This 
aspect of surface texture needs further study. 
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conditions or surfaces exist where the image-intensity equation becomes quite sim¬ 
ple. For example, particular locations on the surface (occluding edge, bright spot, 
hyperbolic region, etc.) may each provide partial information which can be assim¬ 
ilated. Alternately, some common surfaces like water or clouds have very special 
properties leading to quite simple reflectivity functions for certain very common il¬ 
lumination conditions such as extended, hemispheric illumination or direct sunlight 
(see Section 4.2). 
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6.0 Dynamical Properties 

Most surfaces appear in motion on our retina, either because of our own movements 
or because the surface itself is moving. There is a large body of work analysing the 
image correlates of the motion of rigid objects or surfaces (see review by Ullman, 
and Section 5.2 and 5.3). Much less attention has been given to the motion of 
non-rigid objects [although see Ullman (1984), Grzywacz St Hildreth (1987) and 
Koenderink St van Doom (1986)]. Even these studies have not considered the actual 
physics of the deformations underlying surfaces like flags waving in the breeze, or 
fluid motions such as clouds, oil or water. (One exception is Kajiya’s (1984) study 
of cloud formation.) Here we begin an analysis of the dynamics of water-wave 
motion. Then we conclude with a more traditional analysis: the optical flow of 
planar surfaces. Although this study is for a rigid surface, the ambiguity is relevent 
to understanding the appearance of water wave motion, whose image structure 
looks like dark features moving in concert across the image plane. A later study 
of the relation between the smoothness of the (image) velocity field and surface 
rigidity is also relevant, reported here in Section 5.3. This second study shows 
that, in a structural sense, smoothness of the velocity field follows from rigidity of 
motion. Water waves in motion present a smooth velocity field, yet the surface is 
certainly non-rigid. More work is needed to clarify this ambiguity. 

6.1 Non-Linear Wave Interaction Between Long and Short 
Waves (T-.J. Kung) 

1.0 Introduction 

In order to understand the appearance of waves in an image, we explore here the 
origin and behavior of short waves which are less than 40 cm. These waves are 
clearly visible not only to our eye but also to SAR (Synthetic Aperture Radar). 
Surprisingly, SAR images have revealed a wealth of ocean phenomena, some of 
which are correlated to local wind structure, long gravity waves, and large-scale 
circulation patterns. Such findings require a consideration of the interactions and 
energy transfer between wave components having both long and short wavelengths. 

We begin by considering infinitesimal surface water waves, where the differen¬ 
tial equation to be satisfied by the velocity potential is 

*ZS + *MM = 0 (1) 

with boundary conditions 
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. / 2tc 2 * . . 

4>z(x, 0 , t ) = -f— cos — (z - et) 

^t(z, o, t) + o, t) = 0 (*) 

4>z(z, -h, t) =0 

where z-axis is on the undisturbed free surface, z-axis is perpendicular to the undis¬ 
turbed free surface, g is acceleration of gravity, and h is the depth of water. Solving 
the above equations, we obtain the dispersion relations 

w 2 = ffctanh(fcA) (3) 

For the case of deep water, A and tanh(Jfch) « 1. 

ui 2 = gk or c = | = y/gjk (4) 

This is a dispersion wave system. The other limit is that of shallow water, h <SL A 
and tanh(fch) fts kh, therefore 

w 2 = k 2 gh or c = | = sfth (5) 

This is the non-dispersive system. The above analysis is based on the assumption 
that the wave amplitude a is small compared with either depth h or wave length 
A whichever is smaller. Therefore the boundary conditions can be simplified. This 
analysis is actually an approximation of a nonlinear problem in which the nonlinear 
terms are, in some sense, small. Here, we want to examine how the nonlinear terms 
may affect a wave system through the mechanism of wave-wave interaction. 

By way of background, let us consider a very elementary property of oscilla¬ 
tion systems. A linear oscillator subjected to an infinitesimal forcing function is 
described by the equation 

x tt + *?t = ee i0t (6) 

where « is the natural frequency of the system and Q the frequency of the forcing 
function. The response of the system from an initial state of rest is small (of order 
e) unless k 2 is very nearly equal to fl 2 —that is, unless there is resonance between 
the frequency of the forcing function and the natural frequency of the system. If 
the two are precisely equal, the amplitude of the oscillation grows linearly with time 
and becomes arbitrarily large. Indeed, the only way that the response can become 
large is through the phenomenon of resonance. 

A very similar effect is involved in the interaction among wavetrains. Let us 
consider two interaction wavetrains, given to a first approximation by solutions to 
the linear equation. 
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ai exp *'(«!• f-wi*), ^ 

<12 exp *(*2 • * — W2<)» 

where the individual wavenumbers and frequencies are related by the dispersion re¬ 
lation Ui = <•/<(*»)>*’ = 1, 2. The solution to the complete nonlinear equation would 
be expected heuristically to be of the same general form, though the amplitude 
‘a’ may possibly vary with time. We might seek solutions to nonlinear equation 
by substituting expressions of the type (7) into the small nonlinear terms. If the 
lowest-order nonlinearity is quadratic, this substitution leads to expressions of the 
type 

exp t ± * 2 ) - £ - (wj ± w 2 )*] • ( 8 ) 

Now, these terms act as a small amplitude forcing function to the essentially linear 
system and provide an excitation at the wavenumbers (icj ± * 2 ) and frequencies 
(wi ±u» 2 ). The response of the system to this forcing function can be expected to be 
small (that is, of order e) unless resonance occurs—that is, unless the wave r iber 
and frequency at which the forcing is applied correspond to a (wavenumb fre¬ 
quency) pair of a natural wave mode. If they do correspond, we have resonance and 
growth of this new component with continuous energy transfer from the two original 
wave components to this third one. If the initial amplitude of the third component 
is zero, then under resonant conditions it grows linearly until the energy drain 
from the other components begins to reduce their amplitudes and consequently the 
amplitude of the forcing function. 

The above investigation is for discrete waves. A lot of research has studied the 
interaction phenomena between waves having widely different length scales. Since 
the nonlinearity is weak, it takes a long time to see zeroth order changes of the 
system. However, when the group velocity of a wave packet is equal to the velocity 
of a natural mode wave, we expect to see dramatic interchange of energy between 
them just like the case of discrete wave system. Let us consider the infinitesimal 
surface water wave again. If we have a pond with uniform depth h, then a wave with 
wave length A larger than 2 h is considered a long wave which travels at velocity 
y/gK (naturally defined by the system). On the other hand, we can also have a 
wave packet in which the wave lengths are always smaller than h. This wave packet 
is a deep water wave (short-wave) system corresponding to this pond therefore it 
travels at group velocity 

c a = - 7 - = (in linear case) 

v in 2 V k ' 

If the wave number of the short wave is right, the group velocity can be equal to 
y/gK, the long wave velocity. This project is to study the phenomena of resonance 

35 







Richards & Ullman 


Inferences from Images 6.1 


when the group velocity of capillary waves is equal the velocity of long waves. The 
method of multiple scales has been proved to be very effective in studying wave 
interactions and will be used in this analysis. 


2.0 Theory 

For the case of a progressive gravity-capillary wave moving on the free surface of 
a liquid of constant depth h, the undisturbed free surface corresponds to the plane 
z = 0 , where z points vertically upwards, and the bottom is located at z = -h. 
The remaining Cartesian co-ordinates x, y are in the plane of the undisturbed free 
surface, and we choose x to point in the direction of the wave propagation. Since 
the fluid motion is irrotational, a velocity potential <f>(x, y, z, t) satisfying Laplace’s 
equation 

4>xx + <t>yy +<f>zz =0 ~ h < Z < £. (12) 

can be defined, where £(x, y, t) denotes the position of the undulating free surface. 
The boundary conditions for the motion are 

4>z = 0 at 2 = -h. (13) 


4>z = £t + 4>x£x + $yiy at z — £ 


(H) 


g£ + 4>t + - ( <l>x + 4>y + 4>V) 


£xx(l + £ 2 ) + fa y (l + jg] - 2£xy£x£v 

(i + d+tf ) 3/3 


e (15) 


Equation (13) is the boundary condition for the bottom where the normal velocity 
is zero. Equation (14) is the kinematic bound? , -y condition at top surface and 
equation (15) is the dynamic boundary condition. The parameter T is the ratio 
of the surface tension coefficient to the fluid density and g is the gravitational 
acceleration. We suppose that initially (at t = 0) the surface is distorted in the 
manner 


£(x, y, < = 0) = [A(ex, ey)e iKX - A\ <"*} (16) 

where T = n 2 T/g, the asterisk denotes the complex conjugate and « is a nondimen- 
sional parameter measuring the slope of the wavy surface, which has wave length 
2ir/k. The envelope A(c, c) of the surface distortion is allowed to possess a slow 
spatial variation and the frequency is determined uniquely by the value of k through 
the dispersion relation 

u= {?*<T(l+f)}\ (17) 
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where a = tanh(fch). 

We now derive the equations describing the time evolution of A when the 
motion is only weakly nonlinear (0 < e <t; 1). On this basis we assume that (12) 
through (16) have a solution of the form 

* = e *< 1 > + eV 3) +<V 3) + - (18) 

* = ^ 1 >+*V 2) +<V 3) + -- (19) 


Also, we introduce the multiple scales 

£ = «(*- c g t), f) = eg, r = e 2 t, = e 2 (x - c a t ) ( 20 ) 


eg denotes the group velocity and is given by the relation 


c a 


Sw 

6k 


[ a + kk^-a^ f } 

P l 2<r 1 +fJ’ 


u 



( 21 ) 


Substituting these forms into the governing set of equations, solving successively 
the equations resulting from repeated use of the limit process f-»0, with i, y, z, 
t, (, ri and r fixed, and using the notation 

E = exp{t(«cz — wt)} (22) 

one obtains the following results: 

iAr + \A (( +fiA nv = v\A\ 2 A +i/iAQ (23) 


(9h - c 2 )Q it + ghQ vr , = k(\A\\„ (24) 

where 

Q = + {IT? + <, ' (1 “ } |A|2 (25) 

Equation (23) shows that the short wave satisfies the self modulated nonlinear 
Schrodinger equation plus an extra term—the interactions between short waves 
and long waves, while the modulation of long waves is totally excited by short 
waves as can be seen in equation (24). The coefficient v of the cubic nonlinear 
term is singular when Cg(k) = This corresponds to a long-wave/shore- 

wave resonance in which the group velocity of the short (capillary) wave matches 
the phase velocity of the long (gravity) wave. Equation (23) and (24) break down 
under this condition and a different analysis and scaling are required. 

Djordjevic k Redekopp demonstrated that the long-wave/short-wave resonant 
interaction occurs on a much shorter time scale than that for self-modulation. They 
restricted their discussion to the one-dimensional case and introduced the multiple 
scales 
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i = * 3/3 (* - c g t), r = c A > 3 t, fc = c 4/3 (x - c g t),.. . (26) 

The evolution equations describing this interaction turned out to be 

iAr + \A (( = BA (27) 

B r =-o(|A| 2 ) fl (28) 

where 

B((, t) = f*< 0) (29) 

In this project we extend the analysis of Djordjevic Si Redikopp to two-dimen¬ 
sional case. The multiple scales need to be reexamined. It is assumed that <f> ~ 
(1//x) and /x «; 1 since long-wave/short-wave resonant interaction occurs on a 
shorter time scale. Considering equation (23), the second term on left hand side 
and last term on right hand side gives the relationship ~ A<f>( or A ( ~ A<f>. 
Therefore we can write the new streamwise direction scale in terms of old streamwise 
direction scale as {' = {/ft. From first and second term, A r ~ A (( . The new time 
scale is related to the old time scale by r' = r/ft 2 because r' — (£') 2 = (£/m) 3 = 
r/ft 2 . When resonance occurs, c 2 = gh, the first term of equation (24) vanishes. 
The next dominate term is e<f>( T and it has the same order of magnitude as the 
term on right hand side of the same equation. We obtain 

and 

«' = «/,*/> ,,»/»(, 

Equalize the order of second term with the rest of the terms in equation (24), we 
find that the scale for cross-stream direction remains the same, i.e. 

V 1 = *1 = ey- 

As a conclusion, our multiple scales looks like 

i = < 2 ' 3 (x-c„t), f 1 =c 4/3 (x-c a t),... 

T = €*! 3 t 
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Figure 1-4 Numerical solutions for ID equations (27 and 28), showing short 
and long-wave envelopes. 


Based on new scales, we derived the evolution equations describing two-dimen¬ 
sional long-wave/short wave interaction as following: 

iAr + = SA<(> t (32) 

4>t r + Ptvn = -<*(M 2 )* (33) 

If we substitute with 0 as in 1 - D case, i.e. define 0(£, rj, t) = we get 

iAr + AA^ = BA (34) 

Br — 0 J^°BntidZ = —a(|A| 2 )( (35) 

A computer program is developed to solve this set of equations numerically. The 
results are shown in next section. 
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S.O Results 

(A) One-dimensional case. 

Some interesting properties about the 1-D evolution equations [(27) and (28)] were 
found by Koop and Redekopp. The instability is unidirectional in the sense that 
the long wave cannot generate the short wave if the short wave is initially zero. For 
the case of uniform wavetrain with small perturbation, the instability criterion is 
independent of the long-wave amplitude. They also showed that the unstable range 
is 0 < k < 2.182 in terms of the normalized variables where k is the wave number 
of small perturbation. Numerical computation was carried out for two cases (Ref. 
to Figure 1 - 4), k = 2.1 (unstable) and k = 2.25 (stable) with initial conditions: 

5 = 1 + S+e' K * + t-e~ tK( (36) 

L — 0 (37) 

Near-perfect recurrence was observed in unstable case. 

(B) Two-dimensional Case 

The property of unidirection instability, i.e. the long wave cannot generate the 
short wave if the short wave is initially zero, is preserved in 2-D case. However, 
the integral term in equation (35) does make 2-D case behave differently from 1-D 
case. Let us start the discussion with the long term behavior of the long waves. 
From equation (33) 

+ Pinn = -0 (M| 2 )< 

Since the nonlinear effect is weak, we assume that <j> behaves close to linear system, 
i.e. the form of ^ is assumed to be 

* = F(k,l)txp^k£ + ll-jjT^dkdl. (38) 

The stationary phase of 4> gives relationship 

{ = -„ 3 /4r. (39) 

This means that the long waves travel leftward when they propagate outward. It 
suggests that the integration in equation (35) should start from right-most bound¬ 
ary and the integrated value can not be periodical in the stream-wise direction. 
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Instead of using periodical boundary conditions as we did in the 1-D case, we im¬ 
pose zero gradient at the boundaries. Since a uniform wavetrain with periodical 
perturbation can no longer be a solution, we start with a Gaussian as initial con¬ 
dition as shown in Figure (5). As the short wave propagates [Figures (7) and (9)], 
the interaction appears only in stream-wise direction. They never expand in the 
cross-stream direction. This is a proof of unidirection instability. On the other 
hand, the long waves [Figures (8) and (10)] propagate into the left half of the do¬ 
main as we predicted earlier. The boundary affects can be seen in Figure (9) and 
(10) compared with Figures (13) and (14). They both have exactly the same initial 
condition except that the domain of the later case was doubled. 

A couple of different initial conditions [as shown in Figure (15) and (16)] were 
tested. They were so designed that they are close to wavetrain but without the 
periodical property. The results show that the wave number of the perturbation 
does effect the interaction rates, i.e. smaller wave numbers have faster interaction 
rates as in 1-D case. But in all cases, the amplitude of the envelope of short waves 
and the amplitude of long waves only oscillates within a moderate range. The 
unstable case as in Figure (l) was not found. 


4-0 Conclusion 

The equations to describe the two-dimensional resonant interaction between a 
shallow-water gravity wave and a capillary wave have been derived. It is found 
that the resonant interaction occurs on a time scale 0(e -4 / 3 ) just as in the one¬ 
dimensional case. The spatial scale in stream-wise direction is also the same as 
in the 1-D case, namely 0(e -3 ^ 3 ), but the cross-stream direction is different, be¬ 
ing 0(e~ 1 ). The dynamical characteristics of 2-D resonant interaction are quite 
different from the 1-D case, as illustrated by numerical examples. 
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6.2 Optical Flow of Planar Surfaces (S. Ullman) 


I Introduction 

When ulijci l» mow- in 11.. >i.v ... «,ll- i to t lie w«-\wr. tin 

iin:igi-> they ra-u upon I he VM-vw-r - n-iuui undergo ruiuplcx transformations 

tlim-dimensioual (•’•-!)) strmnm- ol tin vi«-\v«-il ol>j*-<is tlu-ir motion in 
space 

A substantia! number of computational studies have investigated this ca¬ 
pacity to recover structure from motion. The main two goals of these studies 
have been (i) to determine the conditions under which the 3-D structure can 
be recovered uniquely from the changing projection, and (ti) to develop meth¬ 
ods for computing the 3-D structure and motion in space from the projected 
transformations. 

With respect to the uniqueness problem, the main conclusion from these 
studies has been that for non-planar rigid objects in motion, the 3-D structure 
and motion are determined uniquely by the projected transformation. Results 
of this type have been established for both perspective and orthographic pro¬ 
jections, and for different forms of input to the recovery process, including 
discrete points in discrete views, discrete points and their velocities, and a 
continuous velocity field. 

The uniqueness results obtained in these studies depended criticially upon 
the object being non-planar. Until recently, the problem of interpreting the 
motion cast by planar surfaces has remained relatively unexplored. The prob¬ 
lem has obvious practical implications, since many surfaces are planar or 
nearly planar, and in these cases methods that assume non-planarity would 
be incorrect or unreliable. As mentioned by Waxman L Ullman (1983) and 
by Longuet-Higgins (1984), the analysis of motion relative to a planar surface 
may be relevant to certain situations of night-landing of an aircraft, where 
the main visual cues are roughly coplanar. 

Hay (1966) was apparently the first to analyze mathematically the visual 
interpretation of moving planes. Hay’s analysis assumed that the visual input 
is given in the form of two discrete views obtained from a set of points in 
motion. 

The main result established by Hay was that the interpretation problem 
in this case exhibits in general a two-fold ambiguity. In addition to the moving 
plane that has actually cause the viewed transformation there is in general 
one additional “confusable" plane. The second solution is entirely different 
from the first in its spatial orientation and motion through space. 

Similar results have been established by Tsai & Huang (1981). based on a 
different method of analysis. They have also shown that when the translattion 
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III 'fi.i' c III iihinp I Ik- normal to ihr tin aml>i>'ml y c!'-.*(>{ wh r- and I hr 

A recent «-i ii< 1\ by Lung'ici-llii-niii' (|>IM) li,i» ulentilii-d an additional 
foiiiliiinii under which die twofold amliinuu di-:i)i)i<ar- lie «t-<> provided a 
different miiln.d iff analysis and a dtllrreui algorithm Tor computing »lit- 3-1) 
parameter*. 

The current paper extend' tin- analysis of ihe planar velocity field (i.e.. 
the velocity field induced by a moving planar surface) in four directions First, 
it uses a different form of input to the recovery process. Instead of discrete 
points, it uses continuous flow parameters such as local vorticity, shear, and 
their derivatives. The use of flow parametsrs in the analysis of the optical 
velocity field was suggested originally by Koenderink L van Doom (1975) 
and by Longuet-Higgins A: Prazdny (1980). An analysis based entirely on this 
form of input has been developed recetly by Waxman Si Ullman (1983). When 
this method was applied to planar surfaces, it was found empirically that 
there wiere in general two distinct 3-D solutions. The result was not proven 
mathematically, however, and the possible existence of additional solutions 
was not ruled out 

An analysis based on input in the form of flow parameters is of interest for 
two reasons. First, as emphasized by Koendrink k van Doom and by Longuet- 
HigginsA: Prazdny. the flow parameters are convenient in the sense that they 
usually have a clear geometric interpretation, and some of them are invariant 
with respect to the choice of a coordinate system. Second, in analogy with 
the non-planar case, it is of interest to examine the interpretation problem 
for different forms of input since the analysis of each case is usually different, 
leading to a different recovery method with somewhat different properties. 

The second direction in which the current paper extends previous inves¬ 
tigations is the consideration of confusable non-planar solutions Previous 
studies have assumed that the moving surface is known to be planar, and ex¬ 
amined the number of possible planar solutions. It is of interest, however, to 
determine the ambiguity of the interpretation when the moving surface is not 
known in advance to be planar. It is shown here that confusable non-planar 
solutions can in general be ruled out. Third, the method of analysis is differ¬ 
ent from previous studies, leading to a different family of possible algorithms. 
Finally, the information available in the orthographic rather than perspec¬ 
tive velocity field of planar surfaces is analyzed. Since under local analysis 
perspective and orthorgraphic projections become almost indistinguishable, 
this analysis indicates some limits on the information that can be extracted 
reliably from a local analysis of the velocity field. 
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2. Tim Wlocity Field of a Moving Surf.irc 


Til.- aual>M-ul tlie planar t.-Wilv li.-U will )....«m tw..M*p* In lln 

step ilii' projected velocity In Id w ill in- cxpre^i-d in term*- of 1 lie .'{-I) shape 
ami million parameter* of tin indin mu siirl.n c. Tin- dev riploin will follow 
the derivation in Wax man <L l liman (I'JSul This step is straightforward, 
and the resulting 2-1) velocity field tsoliviou-ly imiiptcly determined liy thr 
3-D parameters The second siep (Section- 3 1.5) consists of inverting this 
process: given the instantaneous flow Held of a planar surface, the problem is 
to recover the unknown 3-D parameters 

2.1 .Notation 

Following Longuet-Higgins L Prazdny we will use a coordinate system 
(X,Y,Z) moving with the observer relative to the scene. The origin of the 
coordinate system is the vertex of the perspective projection from the scene 
to the image, and the Z-axis is oriented along the instantaneous line of sight. 
The Z-axis intercepts the viewed object at (0,0. Z„). It is assumed that around 
this point the object can be described by a tw ice-differentiable surface (not 
necessarily planar) Z(X,Y). Image coordinates will be denoted by lower-case 
letters (x, y) where x = f, y = \ . V. V, W will denote velocities in space in the 
X,Y,Z directions, and u. v image velocities in the x, y directions. Subscripts 
such as u z ,u y will denote partial derivatives along the x and y direction. 

2.2 The Instantaneous Velocity Field Around the Oriein 

The instantaneous motion of a rigid object can be described by six inde¬ 
pendent parameters. We will denote them by M,,i = 1....6. M\ , Mi, M? are 
the velocities along the X,Y,Z directions scaled by the distance Zo,M t = 
U/Z 0 ,Mi~= V/Z 0 ,M S ■- W/Zo- Mi, M$, Mb are the angular velocities 
around the X,Y,Z axes. The shape parameters we will use are T\,...,Ti, 
the surface orientation at the origin T t - (§j ) 0 , T 2 = (l^) 0 > an< * l ^ e sur ‘ 
face curvature at the origin (scaled b yZo) T s = T 4 = Zo(j?4) 0 , 

Ts = Z 0 (££) 0 . T, are assumed below to be finite. 

The image observables that will be used to recover the unknown 3- 
D shape and motion parameters are the image velocity at the origin (0,0) 
and its first and second derivatives in the i and y directions. Rather than 
using these derivatives themselves, we will use linear combinations of them, 
as suggested originally by Koenderink L van Doom (1975) and by Longuet- 
Higgins At Prazdny (1980). 
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\\r will II'C a srl of-lil iinagr observables denoted b\ (), . O t3 . Tin*) 

art' defined a- follow- 


O. - . O e - 

O, u 3 O.J - r VJ 

O4 - t „ o,„ - t uv 

Or. - J 2 K-r.) 

0 ( . = = O u - u y 


The first two observables evaluated at (0,0) are the image velocity at the 
origin in the z and y directions. The next four can be thought of as describing 
the deformation of a differential neighborhood around the origin. Oj and O 4 
are the rate-of-stretch along the z and y axes, O 5 measures the rate of decrease 
in the angle between line elements oriented along the axes, and 06 is the local 
rate of rotation. O 7 - 0,j are the spatial derivatives of these variables. By 
evaluating explicitly the various derivatives it is straightforward to obtain the 
mathematical relations between the observables 0 ,,...,On and the unknown 
parameters Mi,...,Mo Ti,...Ts (Waxman ii Ullman 1983. See also a similar 
derivation in Longuet-Higgins L Prazdny 1980). The resulting equations are: 

( 1 ) 

O, = -M, - M, Or = -2(M S + M s Ti) + Af,T, 

0 j = — Mi 4 - M4 Og s M4 — M3T2 + Af,T 5 

0$ = M 3 ■+• M\T\ 0o — —Ms — MsT\ + MjT$ 

Oi = M 3 + M,T 2 O,o = 2 (M4 - M,T 2 ) + M 2 T 4 

Os = \(M t Ti + M,Tj) On = |(-M« + MsTs + AfjTj - Af,T 5 ) 

Oo = -A/® + \(MsTi - MjTj) 0 ,j = \(—Ms - MsTi - M1T4 + AfjTs) 

These relations can be used for the recovery of the 3-D shape and motion 
parameters by solving for the M, and T, when the O, are given (measured 
in the image). Longuet-Higgins ii Prazdny (1980) have shown that for non- 
planar patches a similar set of equations has at most three different solutions. 

It was found in computer simulations (Waxman & Ullman 1983) that the 
solution is usually unique, but an analytic uniqueness proof is still lacking. 

We next turn to the planar case and show that eq. ( 1 ) then have in general 
exactly two distinct solutions 

3. The two-fold ambiguity of planar surfaces 

In the planar case the surface curvature parameters Ts,Ti,Ts in ( 1 ) 
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all vanish The resulting equations arc still coupled and non-linear, and the 
liuinlx r of d I'm iI k ! M.iiiliiiii' i- noi Miiiiieiiialeit apparent In l|n- >e<tion n 
is shout) that llirrr are hi general two di<linei solutions. In addition to the 
surface that gave rise to the velocity lield there is in general one (and only 
one) additional surface, engaged in a different motion, that ean produce an 
identical velocity field 


When T 3 = T 4 = 7s = 0 the last four equations in (1) are immediately 
derivable from the preceding eight and can be ignored. If the resulting system 
of eight equations has a solution, it also has a second solution, that is in genera) 
different from the first. This claim is established by giving explicitly a second 
solution in terms of the first. Let (Mi ,..., M 3 ,Ti, 7j) be a solution to (1) (with 
T 3 = T 4 = 7$ = 0). A second solution (Mi, ■ Me,Ti,Ti) can be derived 
explicitly as follows. 

T, = -Mi/M s 
7j = -Mi/Ms 
Mi = -TiMs 

Mi = -TiM 3 (2) 

Ms = Ms 

M 4 = M 4 -Mi- M 3 T, 

Ms = Ms + Mi + M 3 T\ 

Me = M 6 + MiTi - MjTi 

This solution exists provided that M 3 j- 0. Ths case M 3 — 0 is examined 
in section 4 below. The two solutions are dual in the sense that either one 
can be used in (2) to obtain the other. In the two solutions the directions of 
the translation vector and the surface normal are interchanged. If t, n are the 
translation vector and surface normal in the first solution, then the second v 
points in the direction of n, and n in the direction of v. 

We next turn to show that there are no more than two distinct solutions. 
This will be done in two stages. Section 3.2, which is the main step in the 
proof, shows that if two planar surfaces have identical velocity fields then 
they have the same value of M 3 (velocity along the line of sight) Section 3.3 
establishes that there are at most two solutions that share the same value of 
M 3 


3.2 Two planes that induce the same velocity field have the same value of Ms 
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l,i-l * xml A 111-1«o |>laii<-< cngaRcd III nm( lull- ( \/,. .. A/,.) mid (Af ,.Afr.) 

ri-'|n-i iivr !>. '!..ii iiiilmr i.IiiiiuhI n-l..riiy lirliU (i •• identical i»lix rviilil<-< 

III (IJ) Ur i an -—il1 1 it 11.11* 11,< tw,. plain- ii,lrr-r<l Otherwise. T , 

Ty.Tt T 2 . which implies cither (.)(A/,.Af,.) (Af,.Af.,). or (i.) all 

I,f tlx iraii'-laiurn romp..,.ruts Af,. Af;. A/.,. Af,. Af-. A/-,, vanish. This latter 
case, corresponding to pure rotation, is uninteresting since no 3-D informa¬ 
tion is convened by I lie changing image. Lot ( be I lie intersection of * and if. 
This special line in space participates in two different motions and induces the 
same (linear) velocity fields Without loss of generality we can re-orient the 
coordinate system so that the projection of ( coincides with the t axis The 
line f can be described now by Z = ax •* Z 0 . The main step in the proof is 
to consider the intersection line f instead of th>- two planes This line has the 
property that when it participates in motion (Af,, ...,Afc) or (Af j,..., Af 0 ), it 
induces identical velocity fields. 

Consider the situation in which the planes no longer exist, only the single 
line t is moving in space. Let it move with the 3-D motion parameter (Afj - 

Mi,...,Me - Me) From the original coincidence between the velocity fields 

of * and if, and since l lies on both planes, it follows that the velocity field 
projected by ( now vanishes, i.e., it satisfies at the origin the equations: 

v =0 u =0 

ft =0 g. = 0 (3) 

ft =0 

We have transformed the problem of the moving planes into a problem 
concerning a moving straight line. The question is: Under what conditions 
the velocity field of a moving line, as expressed in eq (3), vanishes? In the 
re-oriented coordinate system let us denote 

Af, - Mi — V, Mi - Af, = V v Af s - Af, = V, 

Af« - Af 4 = w x Mi - Me = w y Me - Me = 

The five equations in (3) can be expressed in terms of the six motion pa¬ 
rameters (V,, V f , V,, w y , w t ). The derivation is somewhat lengthy, but straight¬ 
forward. The resulting equations are. 

V z j- Zciu, = 0 

Vy - Z(lW X =0 (A) 

V, - aV, = 0 

aV. + Zrjiiy = 0 

aVy - Zou. = 0 
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From which il follows that: 

l, 0 1, f) 0 l„ 

In term® of the nriginnl plane® t. t. the impliraiion of I . — () in that 
M x — ~M X . In the pammlar roordinate® system in which l project® onto the 
z-axis. it is also true that A/| M \ and A/s - A/?,. 


3.3 The number of distinct solutio n ® cannot exceed two 

Since all the possible planar solutions share the same value of Ms, the 
proof will be completed by showing that for a given A/ 3 there are at most 
two distinct solutions. We proceed along the following plan. If x is a planar 
solution let l be now the intersection line of it with the frontal plane Z = Z 0 . 
We will call t the “tilt line” of *. If we re-orient the coordinate system so that l 
runs along the z-axis, then in the new coordinate system Tj = 0. From eq. (1), 
in this reoriented coordinate system Oj = M3. It can be observed from eq.(l) 
that, if we exclude solutions in which T\ - 0 and M\ - 0 simultaneously, 
the a fixed M 3 and T\ — 0 determine the solution uniquely. W’e will show 
that there are at most two orientations of the coordinate system for which 
O3 — M3. This will imply that (except for two special cases that will be 
examined separately) there are at most two distinct tile lines and therefore at 
most two distinct solutions. We will assume here M s 0, the case Ms = 0 
is examined in section 4. It is convenient for the proof to assume that the 
velocity field satisfies initially O3 = 0 4 . This can be assumed without loss of 
generality, since it is always possible to satisfy the assumption by re-orienting 
the coordinate system (choosing new z, y coordinates) in the following manner. 
Let us rotate the coordinate system along the line of sight by an angle 0 and 
denote the eight observables in the new coordiante system by (Oj,...,Og). 
We will determine an angle 0 such that following the rotation O3 = O4. In 
the rotated coordinate system: 

Os = O s cos 7 0 -f O,sin t 0 L 2OiStn0cos0 
0 4 = O 4 eos 2 0 A- Ossin 7 0 - 2Ossin0cos0 (5) 

to obtain O s = O4, 0 must satisfy 

(Oj - O 4 )cos20 + 2Os6in20 = 0 
assuming 0\ ^ O, 0 is determined by 

fan2/?=2^ (6) 
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There will »lw«\ - |v it solution for i (Not necessarily unique The raw 
O-. - u will lie (\niiiiiieil M'|>aralely) Wc ran h^uiiii ilicreforc llirtl *r have 
mil iei 11\ a riiiiriiiualr in winch <)■■ O, We now roiatc i lie coordinate 

si cm l>\ a new angle o. and denote again I lie observables before the rotation 

lie (O,.0,1 ami following the rotation h> (O,.. .(>.) We v-eh an angle o 

for w hich (>■■ My. In general 0 :t () : rnt~n () 4 stn-n 'JOr.sin aens n 

But since ih. - 0 4 

Of - Of - 20r,sin acos a ( 7 ) 

assuming 0$ * 0 

«tn2a = 

There are four solutions for a in the range |0,2*,, Qj, Qj, Qi *■ x.oj + 
ir. The solution a,,a, + tr are equivalent, they give the same solution for 
(A/,, ...,Mg,Tj,Tj). 

This establishes the claim for the general case Two special cases that 
were excluded from the proof can be analyzed in a similar manner. In the 
first case there is a solution for which T\- 0 and M\- 0 simultaneously. In 
this case the velocity field will have a single direction along which Oj = M s 
Mi, M3, Ms,Ti are then determined uniquely, but there are two solutions for 
M 2 , M 4 , T 2 , in (1) In the second case O s = 0 for every orientation of the 
coordinate system. In this case there are two solutions, in one T t = 7* = 0 
(frontal plane) and in the other M 1 = M3 — 0 (motion along the line of 
sight). In all of these cases the velocity field admits at most two planar 
interpretations. 

4. Unique Solutions 

In the previous section ii was shown that the velocity field of a moving 
plane is compatible, in general, with exactly one additional moving plane. 

There are two special cases under which the two-fold ambiguity disap¬ 
pears, and a degenerate case under which the surface orientation cannot be 
recovered. The discussion of these cases will be brief, since one of these cases 
is discussed in llllman L Waxman (1983) and the other in Longuel-Higgins 
(1984). 

The degenerate case arises when Mj = M j = Ms — 0. The motion in 
this case is pure rotation. The motion parameters can be recovered, but the 
surface orientation remains ambiguous. 

One unambiguous case arises when there is no relative velocity along the 
observer’s line of sight, i.e. M s =0(butA#i and Afj are not both zero) In this 
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ra*«• 1 lie mm imi equal inn- ran l>c -nl\ #*») explii it lx and I lie «o|ijt ion is unique 
'I’ll*' -ecmul Iiiiaml.iiiiimi- . a-.- ..r,-.- vvl.n. M-, ■ u ami l In ..U-erver move- 

.lir.nl> 111\\ .1 rd- ..r ;m.ii .. 11., -urlm .• If lln ..n«ii.al .linn sat.sfie- 

A/, A/ : , Ty . A/j A/:, V-. I lien h ran l.e >«fii from equal imi (2) that I lie 

mtuiiJ -uliiliun ciiiiu ule- with I lit lii-i Tin taller i. m.lil mu i- dimuj-eil in 
Longiiet-Miggiiis. (198-1) lie also suggc-w-d an additional condition that can 
be ti-ied in resolve ilie ambiguity. I>> using an cm ended region of ihe plane 
rallier than the local velocity held. 

In the previous section we have left out the case Os ~ 0. Inspection 
of the equations reveal that this case falls into one of the categories already 
discussed. One case where Os = 0 is obtained when My — Mj — M 3 = 0, 
the ambiguous case discussed above. A second case is when My = Mj = 0 
(but Mi ^ 0) and Ty = Ts = 0. In this case the motion is along the surface 
normal, and the solution is unique. In all other cases there are two distinct 
solutions. 

4.1 Summary 

The ambiguous and non-ambiguous solution can be summarized as fol¬ 
lows. 

1. In the pure rotation case the motion can be recovered but the surface 
orientation remains ambiguous. 

2. If there is no relative velocity along the line of sight {Ms-0 but My M t # 
0), or if the motion satisfied M 1 /M 3 - -Ty, A/j/Mj = -T t (motion 
perpendicular to the plane), then the solution is unique. 

3. In all other cases the local velocity field has a two-fold ambiguity. The 
second solution can be derived in terms of the first by equation (2). 

5. The exclusion of non-planar solutions 

We have seen in the previous section that the velocity field of a moving 
plane has in general one additional planar interpretation. That is, if the 
moving surface is known to be a plane, there are in general two distinct 
solutions The possibility remains, however, that when nothing is known in 
advance about the surface, there are additional non-planar interpretations. 
In this section the question of non-planar ambiguities is considered. In other 
words, the question is whether in addition to the two planar interpretations, 
non-planar solutions are also possible. 

Let *■ be a moving plane and let p denote the vector (M, t Me) of its 
motion parameters. Suppose that the twice-difTerentiable surface P is non- 
planar around the origin and that it induces the same velocity field as ir, and 
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lei )i denote- tin- tiuiinin parameters of P It is shown that if t tie velocity 
field- of /’ and s romi tdc in a neiehliorliood around tin origin, then P is 
in fart planar at tin origin. That is, if T,.T 4 .T : . are the sitrfare oirvature 
paramelers of P at the origin then ]':■ - 7\, - T- = 0. 

W jtliout Jos— of generality we call assume that built surfaces pass through 
the point (11.0. Z..). We ran then dtslinguisli between two cases 

Case J : rr is, tangent to P at (0,0. Z„) In this case it and P have the same 
values for T\.T 2 . From the original equation (1) it can be verified that the 
only ambiguous configuration in this case is when M i and M 2 in (1) both 
vanish. That is. the motion is directed along the line of sight In all other 
cases P must coincide at the origin with *, and satisfy T 3 = Ti — T$ = 0 
Case 2: tr intersects P along some space-curve c. Let -y be the projection of c 
on the image plane Assume first that near the origin 7 is not a straight line 
segments. Longuet-Higgins (1984) has shown that given the image velocities 
of four coplanar points (no three of which are colinear in the image), the 3-D 
motion parameter can be recovered up to the two-fold ambiguity discussed in 
the previous section. 

The implication is that the motion p coincides either with p or with p, 
the motion of the dual solution to In either case we obtain that planar and 
non-planar surfaces with identical motion parameters produce an identical 
velocity field From the original equation (1) it can be verified that the only 
ambiguous configuration in this case arises again when = M 2 = 0. If A/j 
also vanishes, (the pure rotation case), the situation is inherently ambiguous, 
as discussed in the previous section. If M 3 # 0 then * is in fact tangent to P, 
as in the previous case. 

The only remaining possibility is that near the origin 7 is a straight line 
segment. In this case we can use the results of section 3.2. W'ithoul loss of 
generality 7 can be assumed to lie along the x axis. Let M X ,T, denote the 
motion and shape parameters of of P. From 3.2, it follows that 

M 3 = Afj.Afj = Mi,M% — Mi. In addition T t = T\ and 7j = T 3 since 
both are measured along c which is straight line in space common to * and 
P, and 7\ = 0 since w is planar. T 4 ,Ti can now be analyzed using eq (1). 
The equality of the observable O 0 in (1) implies M 2 Ti = M 2 Ti, and since 
Ti = 0, M 2 Ti = 0. Similarly and O l2 together imply M t T 4 = 0, and 
also M\T 4 = 0 since M t = Mj If M t # 0, T 4 = 0 From O !0 together 
with Og it now follows that MjTj = M t Ti and therefore T s = 0. The case 
M 1 = 0 can be analyzed in a similar manner. 

The final conclusion is that T 3 = 7"« = = 0 except for the case 

Mi — M 2 = 0 (motion along the line of sight) 
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In roni liixmii. wr ihii di'tiiigui-li Iwumtii mu cii«e< If tin* relative mo- 
11.>n lirippi-ii* i.i lx- iili.ii* tlx Inn ..f -ifcl.i. 11.< ■- ilu I-xmiI veK.nl> field of a 

live motion parameter;., and whose faugi-nt plane rmririrles with n. Unlike 
the pliiuur two-fold ambiguity this aiitl>is!ijii\ i~ lonil. and can lie resolved l»y 
inspect mg a larger region of the plane 

In the more general rase, in which the motion component parallel to the 
image plane does not vanish, non-plaiiar solutions can he ruled out, and the 
only remaining ambiguity is the two-fold planar ambiguity. 

6. Algorithm 

The proof in section 3 although not entirely constructive, leads to a possi¬ 
ble algorithm for computing the planar solutions. The method will be outlined 
briefly A different algorithm has been suggested by Longuet-Higgins (1984). 

We begin by rotating the coordinate system by an angle 8 to obtain 
Os = O* (where O, are the observables in the rotated coordinated system). 
From (6) we obtain 

tan28 = 

(Provided that Os ^ 0). We can therefore obtain a solution (non unique) 
for ain8,cos8. It is a straightforward computation to then compute the ob¬ 
servables (Oi,...,Og) (Waxman L Ullman, 1983). 

The planar motion equations (first eight equations in (1)) can be viewed 
as eight linear equations in 12 unknown: and X x ,...,Xa, when 

Xj - AfjT,, X 2 = M 2 Tj, X 3 = M t Tj. X 4 = M 2 T,, X s = M*T X , X c = 
M 3 Tj. The eight equations are linearly independent. Furthermore, they can 
be divided into four groups of two equations. 

(M 3 ,X x ,X 2 ) appear only in equations (3,4). 

(Me,Xs,Xi) appear only in (5,6). 

(M t ,Ms,Xs) appear only in (1,7). 

(Mi,M it X a ) in (2,8). 

As a result, each set of unknowns can be solved up to a single scalar: 

(M S ,A,,X 2 ) = (m s ,i,.ar 2 )-t- a(fh s ,5,,x 2 ) 

(Af e ,-Xj,-X«) = (mo.xs,**) + 0(m c ,T3.x 4 ) (8) 

(Mi,M s,*s) = (mi,rns,xs) + "7(^11, ms, is) 

(M 2 ,M 4 ,Xq) = (m 2 , m 4 , x<j) + £(m 2 , m 4 ,x 6 ) 
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i In- hi .in s..j, (i I. ari- «l»*t«Tiimnil l>\ <*ulvit.|> *et>« of Iw.ieq.j*- 

I um* m il.r.-e I...k....u.,- Till- aUr- ... .1. -e n-iiii.ii. i.. I..- di-i (-rifiiiM rt Since 
O, 0 4 . ii f r „,„ (|) i|, M i A/,'/, U 7... i i A, A 2 . ili.Tcfore 

j-| - f>7| = Tj ■ ax* and l.i-nct- a determined. Two possible values for 
an- obtained from the equal ini.-. 

•V, A‘j - A*A« 

A/?.A*, - Afj.Y-. (9) 

A/ 3 A 2 = iVfj.Yo 

respectively. Finally we choose a value for (a,0,~i.6) to satisfy the fourth 
independent relation M t X\ - M iA«. There will be at least one such set 
of values for (a, 0 ,1,6), and at most two. When one solution is known, the 
second can be found immediately using equations (2). 

7. The orthographic velocity field 

Previous sections have established that the parameters of a moving plane 
can be determined up to the two-fold ambiguity from an arbitrarily small 
patch near the origin. When the viewed surface patch is small, perspective 
effects become small, and the recovery process may become unreliable. It 
therefore becomes of interest to analyze the case of orthographic projection 
where perspective effects play no role. (In orthographic projection a space 
point X, Y, Z, projects to an image point z = X, y = V.) 

7.1 The.ort h ogra phic veloc i ty field of no n- planar surfaces 

Let Si be a non-planar surface moving in space. For the orthographic case 
it can be assumed that Sj is fixed at the origin (0,0, Z 0 ) since the translation 
components are immediately recoverable. The rotation of Si can always be 
decomposed into the sum of two components: a rotation with angular velocity 
u- (assumed to be non-zero) about an axis lying somewhere in the frontal 
plane (x-rotation) and a component (z-rotation) about the line of sight Z 
with angular velocity 9. 

A second surface S 2 is said to be a depth tcaling of S t if: 

1 For every point (x,y, z) on Si,(x,y,kz) is a point on Sj for some constant 
k (k* 0). 

2. The rotations Wj and wj are around the same axis, and u>j = w } /k. 

3 . 0 , = 9 2 

For non-planar surfaces the following proposition can be established. If 
Sj is a possible rigid interpretation of a given orthographic velocity field, then 
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V-j i« another possible interpretation if and oil It if it is obtained from .‘>j via 
depth si almp 

Null- that if 0\ II; (' then >1 and |ia\i- dilfirrnl iii-lantaiii-iiii-- 
axes of rotation in spare The urlhupraplm velocity field therefore (!<><•« riot 
determine unique!) the rotation axis m span-. 

Sinre our concern here is primarily with planar surface, the proof will 
be omitted, it can be found in (lllrnan J983) 


7.2 The orthoeraphic velocity field of planar surfaces 

In t he planar case the twofold ambiguity of planar surfaces is combined 
with the inherent depth-scaling ambiguity of orthographic projection As a 
result the orthographic projection of a planar surface admits two interpreta¬ 
tions, each defined up to depth scaling. 

For simplicity of the analysis we can assume that in the planar velocity 
field all the velocity vectors are parallel to the x axis. (If the inducing object’s 
Z-rotation is 8 , then by rotating the observed velocity field by -6 all the 
velocity vectors will become parallel Their direction can be taken as the 
x-axis.) 

The velocity field u(x,y) v(x,y) now has the form: 

«(*, y) = ax + /Sy 

v(x,y)=0 (10) 

If (ui s , w„, «ti ( ) is the angular velocity vector of the rotating surface (as¬ 
sumed to be non-zero) then: 

w v i - w M y = ax + py 

W,X - W x 2 =0 (11) 

One solution to these equations arises when w x = 0. This implies w x = 0 
(if z is not identically zero), and z = ^(w+ 0y). This solution corresponds 
to a plane rotating about the vertical axis. 

If tv. 4 0 then u> x ^ 0 also and z = x. This solution is also a plane, 
with a tilt line along the x axis. 

These two possible interpretations cannot be resolved on the basis of the 
instantaneous velocity field. How much additional information is required 
to guarantee a unique solution? For non-planar objects, it can be shown 
(Ullman 1983) that one additional view is sufficient to remove the depth¬ 
scaling ambiguity For planar objects, the problem is open. 

The orthographic velocity field is thus inherently more ambiguous than 
the perspective one. Instead of two solutions there are two families of solu¬ 
tions, each determined up to depth scaling 


66 









Richards & "Oilman 


Inferences from Images 6.2 


The additional ambiguin of the orthographic veiorit) held implies that 
under local aual>»i-(i.«- twiitg a «iiiall surface patch) the 3-D recovery |>r<Kt-> l > 
i> tint i-iiliri'h -t;il'li- A-pect* of lh« :■.-!) siriirlnrc that are invariant under 
depth scaling are expected to he more stable than others. For planar sur¬ 
faces tliesi no,maiit- include the orii'tilation of the lilt line and the rotation 
romponent around the line of sight Parameters that are not invariant under 
depth scaling such as stirfare slant are expeeled to he less robust. 

8. Summary 

1. The velocity field of a planar surface exhibits in general a two-fold am¬ 
biguity. In addition to the moving plane that has actually induced the 
viewed transformation there is one additional, and in general entirely 
different, planar solution. 

2. There are special cases in which the interpretation of the local velocity- 
field becomes unique. These cases are (i) M 3 = 0 (but M\M 3 ^ 0), 
and (ii) motion directed towards or away from the surface. A degenerate 
case arises for pure rotation. In this case the motion parameters can be 
determined, but the 3-D structure remains undetermined. 

3. Additional non-planar solutions are in general excluded. The exception 
is the case of motion directed parallel to the surface normal. 

4. The two planar solutions can be computed from the eight kinematic ob¬ 
servables, using the algorithm in section 6. 

5 If one of the two planar solutions is known, the dual solution can be 
expressed in terms cf the first using eq 2. 

6 . In the orthographic case there are, instead of two solutions, two families 
of solutions, each determined up to depth-scaling. It is expected that 
for the perspective case only 3-D parameters that are invariant under 
depth-scaling would be robust under local analysis. 

The results explored in this paper are theoretical in nature. They set 
some limits on the performance of any motion perceiving device. It is un¬ 
known, however, to what degree the human visual system can approach these 
theoretical limits It may be of interest, therefore, to test psychophysically 
some of the implications of the above analysis. For example: 

—Can subjects interpret the velocity field of planar surfaces? (i.e., make 
some reliable judgements of relative motion parameters and surface orienta- 
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tion). Re«nli> in (Giboui it. nl 1959) indicate that under <ome condition 
tlu> i' |m»>il>le. altlenit'li tin aii tiracv t* |in.l.i.l)l> mu ver> lne.li 

Can the) intiTpret. at hast to-.. ,|e>>ree. the phmar velocity (ields 

under brief present at ion? If the) do. how do l he) handle the inherent twofold 
ambiguity? 

--('an observers interpret orthographic velocity fields? Can the) recover, 
for example, the tilt of one or both planar solutions? 

Answers to these questions may give us a better insight into the process¬ 
ing employed by the human visual system in the recovery of structure from 
motion. 
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6.3 Rigidity and Smoothness of Motion (S. Ullman and 
A. Yuille) 

The smoothness assumption in measuring visual motion 

The problem of measuring visual motion is, in many situations, undercon¬ 
strained. That is, the information in the changing image is insufficent to 
determine the motion uniquely. This indeterminateness is often referred to as 
the “aperture problem” (for example Marr and Ullman, 1982). An additional 
constraint is therefore required for resolving the ambiguity and determining 
the velocity field uniquely. An important constraint that has been proposed 
for solving the problem is a smoothness constraint. (Fennema and Thomp¬ 
son 1979, Horn and Schunck 1981, Hildreth 1984, Nakayama and Silverman 
1986). When applied to the problem of measuring the velocity of image con¬ 
tours, this smoothness assumption has been formulated by Hildreth in the 
following manner. Of all the velocity fields that are consistent with the trans¬ 
forming image, select the velocity field that minimizes the overall variation. 
That is, if v denotes the velocity at a point, and « the arc length along the 
curve, the preferred velocity field is the one that minimizes the integral: 

/ O’** 

along the curve (the principle of least variation). It has been shown that 
this method of determining the velocity field gives good results under a wide 
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range of conditions, and it seems to correspond in many cases to the velocity 
field perceived by human observers (Hildreth 1986, Nakayama and Silverman 


Justifying the smoothness assumption. 

The main rationale raised in support of the smoothness assumption is that 
the velocity field induced by smooth contours in motion is expected to be 
smooth. This argument is insufficent, however, for justifying the principle of 
least variation. The smoothness assumption implies that if the object is known 
to be smooth, the velocity field should not contain discontinuities. It does not 
imply, however, that “the smoother the better”. Consider, for instance, two 
points in the image separated by a small distance, and moving with image 
velocities Vi and Vj. One may argue that if they lie on the tame surface, 
extremely high values of ||Vi — Vj|| are unlikely, because of some physical 
limitations on the motion of objects in space. But why should one assume 
that a relative velocity of, say, 0.2 deg/sec is less likely than 0.1 deg/sec? Such 
a preference is assumed by the least variation principle (that favors velocity 
fields in which the difference in velocity between neighboring points is as small 
as possible), but cannot be defended only on the basis of a general smoothness 
argument. Clearly, a stronger constraint than just the general smoothness of 
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surfaces is required. 

Another general property that seems to provide a useful constraint in the 
analysis of visual motion is rigidity. Computational studies have shown that 
the 3D stucture of rigid and quasi-rigid objects can be recovered by looking 
for the mo6t rigid interpretation possible of the changing image (Ullman 1979, 
1984). 

In this paper we argue that local rigidity of the object and the principle of 
least variation in the velocity field are related. To investigate this we consider 
a rod moving in three-dimensional space. The rod is allowed to move in a semi¬ 
rigid manner: it can rotate and translate freely and is also allowed a certain 
expansion. We calculate the relative velocity between the two endpoints of 
the rod, as projected on the image plane. The calculations show that under 
a wide range of conditions this velocity distribution is peaked at zero velocity 
and decreases monotonically at higher velocities. As a result, the projected 
velocity field of points linked rigidly together is likely to be consistent with 
the principle of least variation in the velocity field. There are two factors 
contributing to this strong bias towards small differential velocities in the 
projected velocity field: (i) the rigidity of the link between the two points and 
(ii) the effect of projection from 3-D to the 2-D imr • - plane. 

If we imagine an object being made up of a collection of semi-rigid rods 
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this analysis suggests that the projected velocity field is likely to be smooth. 
We conclude that the principle of least variation is a reasonable constraint 
that can be justified for the projected velocity field of rigid or locally rigid 
objects. 

Section 1. The two*dimensional case. 

First we consider a rigid rod moving in two dimensions and being pro¬ 
jected onto a line. Because it is rigid the rod’s velocity can be split into 
rotational and translational components. For our purposes the translational 
component can be ignored as it will not contribute to the derivative of the 
velocity field on the line. 

Let the projected length of the rod be R and the real length be r, the 
angle to the vertical be 9 and the angular velocity be u>. Then the projected 
velocity distribution #(u) is 

♦(ti) = J 6(R- rsm*)*(u - rw cos 9)* r {r)* u (u>)drdudB (1.1) 
where 6 denotes the Dirac delta function, ♦ r (r) and ♦ Ul (u>) are the distribution 
of r and~u/ respectively, raind and rw cosd are the projected lengths and 
velocity respectively. We assume that rods are equally likely to have any 
length between 0 and r m . x and set # r (r) = k between 0 and r maz (with 
kr m „ = 1 ). 
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We can integrate with respect to $ using the first delta function. This 
gives (see Appendix 1) 


$(u) = / - — — - u\/r 2 — R 2 )& u .{u)drdu. (1.2) 

Jr *o Vr 7 - R 7 

We integrate with respect to r to obtain 


♦ ( «) = jT^ 


$ w (u>)du; 
(A 2 u; J + u J )i 


(1.3) 


where aj min = ^!y/r} nax - R 2 - This lower bound comes from the delta func¬ 


tion integration. 

Now we examine the behavior of $(u). If we differentiate it with respect 
to u we get 


*♦ , _ ».(>U) .... 

So, /or any probability distribution $ w (u ;) (which by definition must be a 
positive function) the two terms on the right hand side of (1.4) must be 
negative. $(u) must decrease strictly monotonically with u and has a unique 
maximum at u = 0. 

Thus whatever the distribution of the rotation the meet likely projected 
velocity is zero. This result is similar to that obtained by UUman (1979) 
for the motion of dots. It was shown that if the motion of a random dot is 
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described by an isotropic probability distribution function in three dimensions 
then the probability distribution for the projected two-dimensional motion is 
peaked at 0. 

Note that our result is independent of whether we choose an upper bound 
for r. This can easily be seen by setting r max *-► oo and w 0 in (1.4). In 
fact having an upper bound for r makes dQ/du more negative and makes 
♦(u) decay faster. Henceforth we assume for simplicity that there is no upper 
bound for r. 

To see the connection with the smoothness assumption observe that in 
the limit Hildreth’s smoothness measure can be expressed as 

/»£»’*=E<|>* <«> 

where the sum can be taken over a set of rigid routing rods. The results above 
show that for each rod considered independently the probability distribution 
of the u is peaked at zero. Thus, with this independence assumption, we 
argue that the most likely distribution of the whole contour is the one that 
minimizes (1.5). We will relax the independence assumption in section 4. 

Section 2. The three-dimensional case. 

We now extend the analysis to a rod moving in 3-space. Consider a rod 
of length r projected into the image plane, which has unit normal vector k. 
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The rod’s direction in 3-D space is denoted by r which is a unit vector in the 
direction r. The rod is projected into a vector r p where 


r p = f-(f-k)k. (2.1) 

The rod’s rotation is described by a vector u, where u is the frequency and 

u is the axis of rotation. The velocity v is given by 

r = u x F. (2.2) 

The projected velocity field v v is given by 

v ? = V-(v-k)k. (2.3) 

We consider rods with projected length R and projected velocity 0. We are 
interested in the projected velocity distribution 9(0). Since the projected 
velocity distribution is rotationally symmetric we can express this as U9(U) 
(see Appendix 1). If the rods’ length distribution is $ r (r) and their rotation 
distribution is 9 u (u) then the projected velocity distribution is given by 

*(«) = JS(R- \F p \m - \v,l)* r (f)* u (Z)drd2. (2.4) 

Care is needed in specifying the domain of integration of (2.4). From (2.2) 
we see that only the component of w perpendicular to r contributes to the 
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velocity. We mutt require <3 • f = 0. We impose this by substituting a delta 
function into (2.4) and then integrating over all <3. 

Thus we have 


U*(U) = J 6(R - \? p \MO ■ f)6(U - \* p \)*„{<3)d?du. (2.5) 

We assume all directions of the rod are equally likely and the rods are 
also equally likely to have any length. So we set 


♦(0 = 1 . ( 2 . 6 ) 

We also assume that the rotation is equally likely to be in any direction. Then 
♦ w (<2) is a function of |<2| = u only. 


U*{U) = f 6(R - |r* |)d(<3 • f)6{U - |f,|)* w (w)dfa3. (2.7) 

Choose an orthogonal basis t, J, k. Define angles 0 , <p, a, 0. These angles 
are illustrated in figure 1 . 0 is the angle between the unit vector and the k 
axis. <fi is the angle between the unit vector projected onto the x,y plane and 
the x-axis. Similarly for a and 0. It follows that 

¥ = stnd cotyn + $inB tinip] + eosSk 
<3 = sina co»0i + sina sin0j + eosak 


( 2 . 8 ) 






Figure 1 Dtfiniiiom of the » nglet, tee text. 
Then 

du = du sina dad(3 
df = dr sinO dddp 

and 

'"■* |f^| = rsinB. 

Define the angle between w and r to be r by 













Richards & Ullznan 


Inferences from Images 6.3 


Then from (2.8) we have 


cost = cos$ cos a + sinB sina c os(ip — 0). (2.12) 

Define the angle p by 


cosp — 


(<3*Q g 

luJxfl ** 


(2.13) 


From (2.11) we obtain 


coaptinr as (<3 x f ■ k ). (2.14) 


From (2.8) we find 


cospsinr as *in8 sina(ip - 0) ( 2 . 15 ) 


|tT,| as wsinr fin p. (2.16) 

We use (2.10), (2.11) and (2.16) to write the delta functions as 


*(*-IMW**W*H*M) = S{R~t ain$)6(easr)S{U-Srtintainp). (2.17) 
We now need to change the variables of integration. We have 
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dudf = dudr sina sin$ dad0d8d<p. (2-18) 

The integrand only depends on (p and 0 as <p — 0. We can set V = <P - 0 *od 
integrate out 0. Then we have 


cost = cos8 cosa + sin8 sinacostj) 

cospsinr = sinB sina sinrj) • (2-19) 

dudr = xdudr sina sinB da drj> d8 
Now we must change the variables of integration from o, to r,p: 


/If fe\ 

dadrp = det I I drdp. 

g/ 


From (2.19) we eliminate the terms containing to obtain 


( 2 . 20 ) 


coso = cost cos8 ± sina(sin 2 B - cos V)*- (2-21) 

Define functions F(T,p,a,rp) andG(r,p,a,0) by 

F(r,p,a,il>) = cost - cosBcosa — sinB sinacosij; 

( 2 . 22 ) 

G(r , p ,a, ^) s= cospsinr — sin0sinasinV’ 

It is a standard result of Jacobian transformations (it can be obtained by 

using formula (2.20) twice, once changing variables from a,V> to F,G and 
then changing F,G to r,p) that 
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dadip = det 
From (2.22) we obtain 


(* S?\ (% «\ 

drdp. (2.23) 

V $? «/ vs ? s ?/ 

(9 *\ , 

t I = ain 2 rainp (2.24) 

Vs ? s ?/ 


(% K\ 

det I = stn0s<na{smdcosa - coarpcoa8 aina). (2.25) 

Vs? s?/ 

Thus we calculate 


. j i tinfrtinp . . 

** **nd stna{sindcoso - eoarp coaS sina) iT P ' ( • ) 

We substitute for eos^ from (2.19) and obtain 


. j-a sin'rsinp . . 

***=< Ji7 > 

then using (2.21) we find 


dadrp ■ 


:_ sinrstnp 

*ina(ain 3 8 — co* a p)i 


(2.28) 


Hence 


. .. sinrainpainS , , 

si na atnodadd) = r - r drdp. 

( ain*e-co* 3 p)i 


(2.29) 


ri 








Richards & XJllman 


Inferences from Images 6.3 


We substitute (2.17), (2.18) and (2.29) into (2.7) 


Vi( tt) = J 6(R - rain0)6(eosT)6(u — urstnr sinp) 

- , V ainr sinp sinO . . JA , j 

$ W (w) — -- dr dp dBdu dr 

(stn 2 9 - cos 3 p)i 

Now wetlo the integration with respect to r to obtain 


(2.30) 


U * (V)= I (2.3D 

We integrate with respect to 9 to obtain 


V*(«) - f 


6(U - UT8inp)i u (u)ainp{R/r)dpdudr 
(r 2 - R 2 )b(R 2 /r 2 - cos 2 rtf — 


We integrate with respect to p and obtain 


(2.32) 


R$„(u)<L)dr 


/ r(r 2 - R 2 )k(R?u> 2 + tt 2 - u > 2 r 2 )i(u 2 r 2 - u 2 )* ' ^' 33 ^ 

The domain of integration is restricted to the region of the u - r plane where 
the integrand is real (for simplicity we have omitted the bounds of integration 
in the previous equations). It is specified by 


r>R 

R 2 u 2 + u J > u 2 r 2 


(2.34) 
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Change the variable of integration from r to X where 


wV * u* + X 


(2.35) 


then 


♦(«)=/ 


_ Rv9 w (v)dudX _ 

2XH&u> - *)*(«» + X - + X) 


with domain of integration 


(2.36) 


XKR'J 1 

X>0- (2.37) 

u 3 + X > u 3 R 3 

Define a new variable Y by 

Y 3 = J 3 R 3 -X (2.38) 


then 


♦<«)=j 


_ Ru*J<J)dudY_ _ 

(u*R> - Y 3 )i(* 3 - r J )*(a 2 + h tR 3 - Y 3 ) 


(2.39) 


with 
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Y >0 

Y<vR (2-40) 

Y < u. 

We can further simplify this expression by the substitutions Y = t isinB 
and u = (u/R)Q. This yields 


^ f L (2M> 

It is clear from this expression that 4(u) will decrease with u if and only 
if $u,(u;) does not grow faster than u 2 . Thus for all realistic cases $(u) will 
be a monotonic decreasing function of u. 


Section &. Expansion and Contraction. 

We now consider the situation in two dimensions when the rod ip allowed 
to expand and contract as well as rotate. We can write 

r = rr (XI) 

and differentiating gives 


r = ff + r8b (3.2) 

where n is the unit vector perpendicular to f. We define the expansion coef¬ 
ficient a by 
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f = rj (3 J3) 

and write (3.2) as 

r = rat + r8h. (3.4) 

The projected velocity is then given by 

U = r a aind + rw coad (3.5) 

where 8 is the angle that ? makes with the normal to the image plane and 
8 = u. The projected velocity distribution is given by 

*(U) ss J 8(R - rain8)6(U - (r saind + ruco*8))* t (a)t w (u)drdudsd8 

(3.6) 

where, as in Section 1, we have assumed that ♦(r) » 1 so the rods have no 
preferred lengths. We do the integral of t(U) with respect to 8 to obtain 

*(U)=J(r 3 - R a )- l/7 6(U - s/*-w(r*-* 2 ) 1/2 )*,(s)*„(u;) drdu da. (3.7) 
We now integrate with respect to r to find 

*(U) = J{(U- aR) 1 + u 2 R , )- l ' 2 * t (a)*Mdud3. (3.8) 
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Care most be taken to ensure that this integral is evaluated over the correct 
limits, those for which the integrand is well defined. For u>0we need s < ^ 
and for w < 0 we have s > This domain of integration is shown in figure 
(2). We can split the integral up into two parts 


= + (3-9) 



Figure 2 The iomeun of intefrniion, **e text. 
where 

fumoo ftmU)R 

*i{U)= du «U-nR)* + R i u>)-'''*.(3)* w (")d3 (3.10) 

JvmO Jem- oo 

and 
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* 2 (i7) = r° /**°° a u - * R ) 2 +*v)- , *#.(*)# w («)d». (3.n) 

Jum-oo JtmVfR 
We impose the conditions 

= *.(-•) (3.12) 

and 

#«(«) * #-(-«). (3.13) 

These ensure that contraction is as likely as expansion and that rotation is 
equally likely to be any direction. 

We use (3.12) and (3.13) to rewrite (3.9) as 

rwmoo ftmU/R 

9(U) = J duj {{U- sRf + J2 s w j )- 1/j *.(s)* w (w)d3 

/ w»oo rrm-U/R 

duj ((lf + ^) , + HV)- l / , * > Wi w (tf)(h. (3.14a) 

We now substitute t = As/u and u; = uft/A and obtain 

/"*”<« /'*' ((i - O’ + n’)- 1 ' , *.(»i/«)*„(«n/*)</i 

" Jo-o y t —oo 

+ o ~M + + n 9 r l/ **.MR)KW*)dt. (3.14h) 

Note that if we set 9,(vt/R) = 6(ut/R) we recover (1.3). 
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We differentiate (3.14a) to obtain 


ox rwsoo rt=U/R 

= - du ((U + sR) 2 + R 2 u 2 )~ 3 ' 2 (U + »R)*,(s)* u (u)da 

OU J u -o Jis-oo 

/■w=oo ,t=-U/R 

- du ((U-sR) 2 + R 2 u 2 )- 3 ' 2 (U-sR)* t (s)*„{v)ds. (3.15) 

Jwm 0 J,= -oo 

We can simplify this expression by setting a to a — 2U/R in the second term 
on the right hand side. This yields 


w -- w +wt 

(3.16) 

In the range of integration (ti + sR) is always negative. So the integral will 
be always negative provided $,(s) - $,(s + 2u/R) is negative for s < -u /R. 
If we make s positive, using (3.12), this becomes 


*,(s - 2U/R) - * t (s) > 0, (3.17) 

which is true for all monotonically decreasing functions. Then we have fp- < 0 
everywhere for any function $ w (u»). Hence the projected velocity distribution 
will once more be peaked about U = 0. 

If we relax the conditions (3.17) we are no longer assured of having a 
monotonically decreasing projected velocity distribution for all $ w (u>). If 
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♦,(*) is linear in a we obtain a flat distribution 

*(U) = const. (3.18) 

If 4«(s) is quadratic in a we find 

= -2 J + jj *«!«*(.*+w>)-V> !£#„<«) (3.19) 

which will only be negative for some distributions $u,(u;). 

Thus allowing the rod to expand and contract reduces the strength of the 
result of Section 1. But if we make the reasonable assumption that *,(s) is 
at most linear in s the result will still hold. 

Section 4. Curves from rods. 

The situation becomes more complex when we consider several rods 
joined together to form a curve. The results above show that for each rod 
individually the projected velocity distribution is peaked at zero. If the rods 
are joined, however, these probability distributions are no longer independent. 
We now sketch an argument suggesting that this loss of independence does 
not affect the result. 

Suppose we have a set of probability distributions Pi(x<), « = l,..N 

which are all peaked at zero. We now define the joint probability distribution 
P{2) where 2 = (*j,xj,...,**)• If we impose no constraints this distribution 
is given by 
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P(x)dx = n i P i (x i )ll i dx i . (4.1) 

Now impose a consistency constraint x, = 0. We now have 

j P(x)dx = J IliPiixiMY,**)* 1 ** 1 **' (4.2) 

The Xi’s can no longer vary independently and the delta function imposes the 
consistency constraint. We can get integrate out the delta function thereby 
reducing the number of independent variable to N — 1. It is convenient to 
choose these to be x< - x i+ i, for i = 1,..., JV — 1. We can now write (4.2) as 

J P{2)dZ = J n iPiixi - x i+1 )n id( Xi - Xi+i). (4.3) 

We can set yi — x, - x,+j, » = 1, ...,N — 1 and rewrite 

P(if) = Pi(yi)P2(y7)-Ps-i(yN-i)Ps{-yi - ya -... - vs-i)- (4.4) 
It is straightforward to see that the maximum value of P(y) occurs when the 
Vi are all zero. Thus P is still peaked at zero even with the constraints. 

Section 5. An alternative approach. 

The results derived in the last few sections support the assumption (Hil¬ 
dreth 1984) that the velocity field along a contour is smooth. This assumption 
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is used to solve the aperture problem by prescribing the smoothest possible 
velocity field consistent with the data. 

An alternative approach to the aperture problem was taken by Waxman 
and Wohn (1986). They showed that if the object was locally planar and 
moving rigidly then, assuming perspective projection, the local velocity field 
on the image plane would be locally quadratic in the x,y coordinates. They 
then defined regions in which these expansions were valid and checked for 
consistency between these regions. In these regions they could do a least 
squares best fit to find the coefficients of these quadratic polonomials, and 
hence the velocity field. 

This approach, however, can also be used to justify the smoothest velocity 
field assumption. The orthographic projection of the motion of a rigid body 
will obey 


Vg = A + By + Cz, 


(5.1a) 


v v = D + Ex + Fz , (5..*6) 

where A, B,C, D, E, F are constants and z is the (unknown) depth coordinate. 
If we also assume the object is planar then z obeys 
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z = px + qy + a, (5.2) 

where p,q t r axe constants. Substituting (5.2) into (5.1) gives linear expres¬ 
sions for Vg_ and in terms of x and y. So if we assume local planarity and 
local rigidity (in the style of Waxman and Wohn (1986)) the velocity fields 
will locally be lm«*r in x and y. The measure of smoothest velocity, /(v), 
used by Hildreth is the integral along the curve of the function 




dv dv 
da da' 


(5.3) 


where s is the arc length. The velocity fields will be locally linear if and only 
if J(v) is zero. Thus we can think of the smoothest velocity field approach as 
a local method of assuming local rigidity and local planarity. 


Sectio i 0. Conclusion. 

The arguments in the first four sections of this paper suggest that locally 
rigid, or semi-rigid objects, will tend to project a smooth velocity field in 
the image. Moreover, assuming random motions and limited expansion or 
contraction, this field will tend to be a. smooth as possible. In the fifth 
section we noted that rigidity of an object and local smoothness of its surface 


will also lead to a smooth image motion. 
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Tnese arguments support the view that maximixing smoothness is a good 
heuristic to use for motion correspondence and that it is a sensible way to solve 
the aperture problem. 
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Appendix 

We first describe the general method for integrating delta functions. 

Suppose we have an integral I(a,b) 

I(a,b) = J*f(z)6(z-z 0 )dz (A. 1 ) 

where 6(z) is the Dirac delta function, /(z) is an arbitrary function, zo an 
arbitrary point and b > a. The value of the integral is 

J(°,6) = /(z 0 ), if z 0 c[o,6] ( A.2a ) 

/(o,6) = 0, otherwise. (A.26) 

This result can be generalized to integrals of form 
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J(o,6) = J* f(z)S(g(z) — c)dx {A. 3) 

where g(x) is an arbitrary function and c an arbitrary number. All points Xi 
with g(xi) = e will contribute to this interval. SnppoBe there is only one such 
point. If there are several we can divide the integral up into regions with only 
one such point. Consider one such point * = 0. The function g(x) can be 
expanded in a Taylor series about this point 

<K*) = c + s'(0)(*) + O(x 2 ). (A. 4) 

If we change the coordinate to u where u = g*( 0)x we can write the 
integral as 


/(a,6)= /'‘%«/s’(0))«(« + O(« ! ))-£L. (.4.5) 

•//(Oja ^(0) 

The value of the integral will depend on the sign of tf'(O). If it is negative 
the bounds of the integral will be reversed and the ! ' <^ral will change sign. 
Therefore 


■' < ° , ‘ ) = /( 0 ) isW’ 0<M1 {Aia) 

J(a,b) = 0, otherwise. (A.6b) 
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We now consider the form of the probability distribution function for a 
rotational symmetric vector in two-dimensional space. Suppose the function 
is ia{Q)du. If it is rotationally symmetric (and hence depends only on the 
modulus u of u) it can be written 


$2(Q)du = $ u (u)ududp, (A.7) 

by changing to radial coordinates u,<p in the u> space. The (p component can 
be integrated out to give a distribution 
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7.0 Assessment and Future Directions 

Computer Graphics provides a powerful too! for understanding vision. Although it 
is exceedingly time consuming to generate realistic models of a scene from physical 
principles , the rewards are great. There is no other method available at present 
of comparable power that allows us to isolate the physical parameters relevant 
to “seeing” natural surfaces or objects. Renditions now exist for many classes of 
plants and trees, and for several kinds of natural phenomena like water, clouds and 
sunsets (see SIGGRAPH and Fournier ii Reeves, 1987). These programs present 
an unusual opportunity for the visual psychophysicist who wishes to understand 
the perception of natural surfaces, and how inferences about material types and 
objects can be made from images. 

One way to further the Graphics Psychophysics approach would be to encour¬ 
age psychophysicists to collaborate with those laboratories actively engaged in the 
generation of natural phenomena using graphics. Typically these graphics labs are 
interested in selecting the best rendition parameters, but do not bother to conduct 
parametric studies which the psychophysicist is trained to do. Graphics laborato¬ 
ries at Cornell (Architecture), University of North Carolina, Toronto, Cal. Tech., 
NYIT and MIT-Media might be approachable hosts. A special sabbatical or post¬ 
doc fellowship to a young investigator would encourage this collaboration at little 
cost. Psychophysical studies of the rendition of even a half-dozen different materials 
would greatly advance image understanding. 

Finally, this project only began to explore time-dependent phenomena, namely 
the motion and creation of water waves, plus some semi-rigid motions. This rel¬ 
atively unexplored area of non-rigid motions should be contrasted with the large 
number of studies which are devoted to rigid body motions or to flow fields. Non- 
rigid motions include not only the wave motions of fluids such as water, honey, 
milk, oil, etc. (which we can tell apart psychophysically), but also turbulent mo¬ 
tions, deformations, and articulated motions. Here we also have lacunae well worth 
filling. 
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Appendix I: Relating Fractal Surfaces to Fractal Images 

Our aim is to recover a 3D fractal measure of surface roughness. As a first step, we 
need to show that the resultant image contains a measure related to the 3D fractal 
specification for the surface. Pentland (1983) and Kube k Pentland (1986) have 
already shown that a convex Lambertian 3D surface of constant albedo with a spa¬ 
tially isotropic fractal Brownian shape seen under constant illumination produces 
an image whose intensity surface is also fractal Brownian with a fractal dimension 
identical to that of the components of the surface normal. Specifically, given 

/(*,») = pE(NL) (1) 

where p , E, L are constant and IV is a Brownian fractal, then /(x,y) will be fractal. 

We wish to prove a similar result for glossy surfaces. In particular, we require 
that I will obey the rule 

V / 

where h is the Hausdorf dimension and 

I(x,y) = E(NH)F fieaael (3) 

where H = L + V (i.e. bisector of illuminant and viewer directions). F freimel is the 
Fresnel reflectance function. Unfortunately, F freanei is dependent upon the relation 
between N and L (incident angle), even when N and H align, which is seldom for 
point sources. Together, these factors generally will prohibit a point source from 
generating specularity regions that satisfy a fractal image statistic. 

However, most scenes are illuminated also by diffuse, hemispheric sources. 
Under these conditions, disregarding occlusions, for every microfacet viewed there 
will be a patch of the hemisphere reflected off the facet. Equation (3) now becomes 

I(z, V ) = E*F ft " Del (NV) (4) 

where the Fresnel term is a function only of the emittant angle, which equals the 
incident angle, and E* is illumination per solid angle. 

We next approximate the Fresnel reflectance by 

F = (1 — V-N) n = (1 — cos 9) n (5) 

where n is a parameter related to the index of refraction, and controls the rate at 
which the Fresnel reflectance falls toward zero as the light beam moves away from 
its maximum reflectance of 1.0 for the grazing condition. For incident angles 6 
greater than 60° we can approximate (5) by its Taylor Series expansion, retaining 
only the first terms: 
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But as long as V is a constant, which it is for viewing glossy surfaces at sufficient 
distances, then the fractal nature of N will result in a scaled version of (8) satisfying 
condition (1). Hence glossy reflections off a fractal 3D surface seen under hemifield 
illumination will produce a fractal image-intensity surface. 
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Appendix II Recovering Material Properties from Sound 
(R.P. Wildes) 

1.0 Introduction 

Physical modeling of world events occupies an important position in. computational 
approaches to perception and motor control (Brady et al., 1982; Horn, 1970; Marr, 
1982; Richards, 1988). Such models are dual in purpose. First, they lend to a precise 
statement of the problem under consideration. Second, they expose constraints on 
the solution space, which indicate how the problem can be solved. Here, we use this 
approach for an initial attack on a previously unsolved problem in audition: How it 
is possible to recover the material property of an object from the sound generated 
when it is struck? 

Consider some examples. When you hit a drinking glass with your finger-nail 
it produces a distinct “ringing” sound. In comparison, if you strike a log with your 
knuckles it gives off a short “thud”. In either case you easily identify the material 
of the struck object. There are other cases where the sound from striking an object 
leaves one with no clear impression of material type. Objects which are clamped or 
otherwise artificially damped often belong to this latter class. Our goal is then two¬ 
fold: First, we seek to discover a physical parameter of the sound following impact 
which is intrinsically related to material type. Second, we desire to recognize those 
situations where such material identification is not possible. Guided by a model 
of vibrating solids we shall choose to satisfy these goals with a measure of energy 
dissipation during vibration—an intrinsic material property. 

2.0 A Physical Model 

When a solid object is subject to an impact, much of the resultant sound is due to 
the vibration of the object. Therefore, we study the mechanics of vibrating solids, 
and take that motion as analogous to the sound production. In the following 
presentation we do not consider the initial transient due to impact. Instead, we 
concentrate on the steady-state and damped behavior of the solid. We proceed in 
two stages. First, we introduce some basic concepts concerning the deformation of 
solids. Second, we develop a particular model of a vibrating solid: the standard 
anelastic linear solid described by Zener (1948). 1 

1 Many of the results in Sections 2.1, 2.2 and 3.1 are known in the literature on anelas- 
ticity. For further details see Nowick it Berry, 1972; Wert, 1986; Zener, 1948. 
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2.1 Basic Concepts 


Hooke’s law states that an ideal elastic material can be described by the relation 
o = Me (1) 

Where: a is a stress variable (corresponding to a force); e is a strain variable 
(corresponding to a displacement); and M is an appropriate modulus of elasticity. 2 
In many situations it is convenient to define J = -fa, ns the modulus of compliance. 
Then (1) becomes 

t - Jo (2) 

Elastic behavior so defined is seen to be characterized by three conditions. First, 
there is a 1 : 1 correspondence between stress and strain. Second, the stress-strain 
relationship is linear. Third, there is no time dependence in the mapping between 
stress and strain. Clearly, any real solid is only approximately characterized by 
these conditions. In most cases, condition three is a particularly erroneous assump¬ 
tion. Therefore, we proceed by considering cases where there is a time dependence 
between stress and strain. Solid materials characterized by such a time dependence 
are called anelastic solids. 

Suppose a specimen is subjected to a periodic stress 

<t = (3) 


Here: a o is the stress amplitude; u = 2xf is the angular frequency; t is time; and 
« = Then, our assumptions of linearity and time dependence dictate that 

strain has the form 


(<) 


with eo the strain amplitude and <f> a phase angle. Then, 4> captures the notion 
of a time dependence between stress and strain. For subsequent operations it is 
convenient to separate (4) into its real and imaginary parts. Letting be the 
component of strain in phase with stress and e 2 be that component of strain § out 
of phase with stress we have 

« - ( £ l ~ «2)e‘"‘ (5) 

Now, in analogy with (2) we define the complex compliance 

J(w) = ^ (6) 

Then, dividing (4) by (3) gives 

2 Notice that for a fine grained analysis o and e would be second order tensors. For our 
purposes it suffices to consider them as scalars. 
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= (V 

<To 

The complex compliance can also be represented as a sum of real and imaginary 
parts by dividing (5) by (3) and defining 

•J(w) = Moj) - ij 2 (u) (=) 

We proceed by expanding (7) in terms of Euler’s relation 3 and then equating 
read and imaginary parts of (7) and (8) to get 

Ji(u) = — cos 4>(u) 

2 <» 

<tq 

From which we reach the important result that 

tan * = ( 10 ) 

In the anelasticity literature tan <f> is referred to as internal friction. Internal friction 
is an intrinsic property of a given material; it measures the degree of anelasticity. 
We shall ultimately recover a measure closely related to internal friction as our 
characterization of solid materials. 

Finally, it is worth noting that analogous calculations can be carried out in 
terms of the modulus of elasticity. Such derivations would begin by defining 

*M = f (ii) 

with the ultimate conclusion that 

tan ^ = XT ( 12 ) 


2.2 The Standard Anelastic Linear Solid 

The model of solids that we wish to consider is depicted as a mechanical structure 
in Figure 1. We refer to this model as the standard anelastic linear solid. The 
model consists of three elements. Elements a and b are Hookean springs; let them 
have corresponding compliances J a and Jj>. Element c is a Newtonian dashpot, 
which serves to provide internal friction. A Newtonian dashpot can be thought of 
as a plunger moving in a viscous liquid. Its velocity of motion is proportional to 
the applied force. In terms of stress and strain 

<r = i»< (13) 

3 e' 8 = cos 6 + * sin 6 
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Figure I A mechanical model of the standard anelastic linear solid. J a and 
J b are compliances of the two springs, n is the viscosity of the “liquid* in the 
dashpot 


with T) the viscosity of the liquid. For our purposes we let n = , where t„ is 

the time required for the system to achieve equilibrium when a constant stress 
is applied. Notice that this model captures our intuitive notion of what happens 
when a force is applied to a solid: Initially the solid deforms (element a). With 
time further deformation takes place (element c). When the force is removed the 
solid returns to normal (elements a and b).* 

In order to be more precise, we characterize the behavior of our model in terms 
of the differential equation 

ao<r + ai& ~ boe + bii (14) 

where the derivatives are with respect to time. We begin by setting 


c a = Ja^a 

fb = Jb^b ( 15 ) 



4 In some cases the solid never returns to its original configuration. This introduces the 
concept of visco-elasticity; we shall not consider such cases. 
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Figure 2 A second mechanical model for the standard anelastic linear solid. 


We now introduce the standard rules for combining multiple stresses and strains: 
For series combination, stresses are equal while strains add. The reverse is true for 
parallel combinations. Applying these rules to the model of Figure 1 gives 


c = e a + *6 
ffc = Cc 
O = <7 a 

= ± (16) 

= <T b + O e 
_ Cfe + r 9 k c 

Jb 

We can rearrange (16) to the form (14) as 

J<r + T 0 J a & = € + re (17) 


where J = J a -+■ Jb. 

We now characterize our model in terms of the standard measure of internal 
friction tan<£. Following methods similar to those of Section 2.1, we substitute (S) 
and (5) into (17) and equate real and imaginary parts to find that 


1- (wr «) 2 


H = Jb~ 


■ (u>To ) 2 


(18) 


Then, upon substituting into (10) we get 


tan <f> = 


Jb 


J + J.(wnr) 2 


(19) 
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We continue by deriving a more convenient form for tan <f>. Consider the model 
of Figure 2. This model can be shown to be formally equivalent to the model of 
Figure 1. It is customary to represent this latter model in terms of the elastic 
moduli and the time to strain equilibrium r e . Following the methods used above 
and using analogous symbols, we can find that 




Ma + M b (uT t y 


M-2 — Mb'. 


f- (wn ) 2 

<JT t 


( 20 ) 


(- (wr e ) 2 

Finally, to get our simpler expression for tan^ we combine (18), (19) and (20) to arrive 


tan <j> 


J b _ u/r 

(, JaJ ) 2 1 + («r) 2 


( 21 ) 


with r the geometric mean of r„ and r e . 

We conclude this section with a recapitulation. In this section we began by 
introducing the concepts of stress, strain and internal friction. We proceeded to 
present the standard linear model of anelastic solids in terms of these three param¬ 
eters. In Section 3 we make use of these relations to classify solid materials from 
the sound they produce following impact. 


S.O Recovery of Material Type 

When an anelastic solid is struck and set into vibration its behavior, and hence the 
sound it generates following a brief transient, is dictated by our standard linear 
model. Our goal is to extract some intrinsic parameter of the solid material from 
this dynamic behavior. As noted earlier, internal friction, tan <f>, is just such a 
measure. We now proceed to derive two measures of internal friction and then 
show how they can work together. 


S.l Internal Friction and Bandwidth, Q~ 1 

We begin by relating internal friction to peak vibration of our standard linear 
model. To facilitate this set of derivations, we follow (Zener, 1948) and momentar¬ 
ily consider the mechanical system as having only one degree-of-freedom. This is 
depicted in Figure 3. Here we have lumped the inertia of the system into a single 
member I, while the anelastic and elastic components are lumped into an anelastic 
spring with complex modulus 6 

6 It can be shown that an arbitrary standard anelastic linear system in vibration can 
always be so reduced, see (Nowick k Berry, 1972, appendix A.) 
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Figure S A lumped model of Figures 1 and 2. 

M = Af(l + iUn<f>) (22) 

From Newton’s law we observe that 

n = <t ( 23 ) 

Assuming periodic motion, we also have that 


and as always 

<r — Me 

We can solve for e by substituting (22), (24), and (25) into (23) to get 

e — -r- 

(i ~ i ir) + * tan ^ 

By inspection, we find (26) at a maximum when 



(24) 

(25) 


(26) 


(27) 


Further consideration shows that (26) is at 2~ $ maximum amplitude when 

^ = 4a = y(l±tan^) (28) 

We arrive at our first measure of internal friction by noting that 

-*■ = Q := tan 0 (29) 
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where we use Q in accordance with the corresponding concept from electrical circuit 
theory. Thus, internal friction can be measured in terms of the sharpness of the 
peak around the maximum of vibration of our specimen. This will correspond to a 
sharpness in the peak of the acoustic signal generated by this vibration. It is this 
acoustic peak which we suggest measuring. 


S.S Internal Friction and Decay Rate t e 


We now seek a measure of internal friction in terms of the decay in amplitude of a 
vibrating anelastic solid. Combining (22), (23) and (25) we find that 

It + M(1 +1 tan 4>)t = 0 (30) 

describes the behavior of our standard linear system as the vibration decays. The 
solution to (30) is found to be 

< = eoe~ Aft e iut (31) 


with A — ln(A n fA n +i), where A n and An+i are successive peaks in the amplitude 
of the decaying tone of frequency /. (A is called the logarithmic decrement.) We 
solve for the internal friction by substituting (31) into (30) and equating imaginary 
parts to find that 

A 


= tan^ 


(32) 


Because of the potential difficulties of measuring two successive amplitudes, we 
proceed a step further. Letting t e denote the time required for the amplitude to 
decrease to £ of its original value, we note that 

(33) 


e / A wA 

Rearranging and substituting (33) into (32) then yields 

1 2 
tan 4> = —rr — ~T 
nfU ut e 


(34) 


Our second measure of internal friction is then given in terms of the time it takes 
the amplitude of vibration (and therefore of the sound) to decrease to some fraction 
of the original amplitude. 
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S.S Specific Loss 

We now have in hand two measures of internal friction (29) and (34). Examination 
of these relations reveals that tan ^ is a function of frequency. For our purposes, 
it would be more convenient to have a single characterization of anelastic behavior 
rather than an entire “damping spectrum”. Such a measure is available in terms of 
loss per cycle, known as “specific loss”. To motivate this measure we notice from 
(21) that for (ur) 2 > l, wtan <j> is essentially independent of frequency. Empirical 
investigations of solids (Bennewitz k Rotger, 1936; Gemant k Jackson, 1937) con¬ 
firm that this relation is true for a very wide range of vibrational frequencies and 
materials. Therefore, we define specific loss as 

f = wtan^ (35) 

as our measure of anelasticity and as our key to classifying solid materials by sound. 


S-4 A Constraint Curve 

The final task we face is to develop a method of telling when our two measures 
of internal friction will yield an accurate characterization of a given sample. For 
example, let us say an object is struck and our auditory syBtem recognizes its pitch 
as Ag = 880 Hz with a decay rate of 1 sec for t e . By equation (34) we can recover 
a measure of internal friction, and hence characterize the material. But how do we 
know this answer is correct? To confirm that the material behaves as our standard 
anelastic linear solid, we need a second, corroborating measure. But of course, this 
can be obtained by recovering the bandwidth, Q" 1 , of the sound about 880 Hz, 
using equation (29). In this case, if each measure does not yield the same value for 
tan 4>, then we must disregard this characterization. 

By clotting t e versus Qu -1 as shown in Figure 4, we can construct a constraint 
line tha. aolds for materials obeying our standard anelastic linear model. The line 
is independent of frequency because each measure of internal friction has been 
recast as “specific loss” (35). Clearly, all valid measurements must lie along this 
constraint curve. Here we show how particular material types spread out along the 
curve with metal at the upper right and rubber at the lower left. The data have been 
collected from Gemant k Jackson (1937), Berg k Stork (1982) and the Handbook 
of Chemistry and Physics. Exactly how fine a distinction can be realistically made 
along the constraint curve is largely an empirical question (Waller, 1938; Warren 
k Verbrugge, 1984). 
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Flgnre 4 The constant curve relating the two measures of anelasticity for our 
standard linear solid model. Decay: t e ; tuning width: Q~ l ■ 


4-0 Summary 

We have described how to make auditory measures of an intrinsic parameter of solid 
materials. The parameter chosen is called specific loss, and is related to internal 
friction, a measure of anelasticity for a given material. The two measures presented 
are related to the width of the resonant peak of a sound and the time of its decay. 
Crucial to our development has been an understanding of a standard anelastic linear 
system, a physical model of the mechanics of solid vibration. We hope this work 
demonstrates the usefulness of using an explicit physical model of sound production 
in attempts to study sound recognition. 
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